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conducted by the Electrical Engineering Department of Auburn University 
during the period 17 June 1965 through 15 October 1972. This study 
report completes Contract No. NAS8-20163, granted to Engineering 
Experiment Station, Auburn, Alabama, by George C. Marshall Space Flight 
Center, National Aeronautics and Space Administration, Huntsville, 
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SURVEY OF DIGITAL FILTERING 
H. Troy Nagle, Jr. 

ABSTRACT 

A three part survey is made of the state-of-the-art in digital 
filtering. Part one presents background material including sampled- 
data transformations and the discrete Fourier transform. Part two, 
digital filter theory, gives an in-depth coverage of filter categories, 
transfer function synthesis, quantization and other non-linear errors, 
filter structures and computer aided design. Part three presents 
hardware mechanization techniques. Implementation by general-purpose, 
mini-, and special-purpose computer are presented. 
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I. INTRODUCTION TO DIGITAL FILTERING 

Digital Filtering may be described as the process by which input 
discrete-time sequences of numbers with discrete amplitudes are trans- 
formed into output discrete-time sequences of numbers with discrete 
amplitudes. The transformation process (the digital filter) may be 
described as a set of difference equations which may be progra mme d on 
a general-purpose computer, or realized with specially designed devices. 

The Computer Model 

An example digital filter is shown in Fig. 1. The input signal 
e i (t) is in analog form and is sampled every T seconds by an Analog- 
to-Digital Converter (A/D). The input samples e^nT), n an integer, are 
in binary 2’s complement form and are supplied to the computing device, 
which may be a general-purpose or special-purpose computer. The com- 
puting device is programmed to calculate the filter output samples e 0 (nT) 
which are fed to an output hold register. This register may actually be 
considered to be part of the computing device. The output samples e 0 (nT) 
are held in the register until a new output sample is calculated and sup- 
plied to the output hold register. The D/A converter produces an ana- 
log output signal e Q (t) whose characteristic form is shown in Fig. 2. 

The digital filter of Fig. 1 operates as follows: A pulse from the 

digital filter control unit at t=nT instructs the A/D to calculate 
e^(nT). However this sampled value of e^(t) is not available to the 
computing device until time nT + T a, where Ta is the total A/D 
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conversion time. Once the input sample has arrived at the computing 
device, the output sample is calculated and sent to the output hold 

device at t=nT+TafTc, where Tc is the computing time. Thus the output 
e 0 (nT) is actually e 0 (nT+Ta+Tc) . With modem technology, Ta+Tc can be 
designed to take less than lys. Hence for sampling rates of up to 
200KHz(T=5ys) , T is much greater than Ta+Tc and hence, 
e 0 (nT+Ta+Tc) ~eo( n T)» 

The z -Domain Model 

The digital filter of Fig. 1 has been examined and described from 
a hardware or functional point of view. A mathematical model for this 
filter is shown in Fig. 3 which employs the well known z-transform. 

Fig. 3a demonstrates a discrete time model for the digital filter where 
the switches labeled T represent impulse samplers and the block labeled 
G^ 0 (s) represents a "zero-order hold" device. Fig. 3b illustrates the 
transfer function representation of the computing device itself. Fig. 

3 differs from Fig. 1 in that the computing device of Fig. 1 uses the 
amplitude of the input (and output) samples to calculate new output 
samples e Q (nt) ; Fig. 3 uses impulse functions weighted by the amplitude 
of the input (and output) samples to calculate new output impulses 
e Q *(t). The impulse samples, zero-hold, and z-transform will be dis- 
cussed in more detail later. 

Scone of Digital Filtering 

The digital filter models presented above were for conventional one- 
dimensional processing of a single input variable. Although this con- 
cept of digital filtering is most widely accepted, many other researchers 
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have applied the label to more general schemes. Digital filtering in the 
last few years has come to also mean optimal state estimation, discrete 
Fourier transformation, high speed convolution, non-linear discrete 
filtering, two-dimension image processing, random and multirate sampling, 
block recursion, least-mean squares filtering, quantization optimization, 
computer programming, and hardware implementation. All of these topics, 
and others, will be introduced in what follows. Emphasis will be, how- 
ever, on the standard case of linear, one-dimensional digital filtering. 

A prerequisite to understanding the theory of digital filtering is a 
mathematical background in sampled-data transforms, Fourier transforms, 
convolution, discrete state variables and stocastic processes. Hence, 
these topics are now reviewed briefly. 



II. STANDARD Z-TRANSFORM 

The most common sampled - data transformation is called the standard 
z-transform. It is used to describe both the sampling process for the 
digital filter input signal and the discrete transfer function of the 
filter itself. 

Impulse Sampling 

The z-transform is used to represent mathematically a discrete- 
time system. The discrete-time intervals are produced by periodic 
impulse samplers. Consider Fig. 4. The Laplace transfer function G(s) 
is an analog filter; its input is a Dirac delta function. The filter 
output g(t) is periodically impulse sampled every T seconds. The 
sampled output may be expressed as 


g*(t) = g(t) E 6 (t-kT) 

k=0 


= E 
k=0 


g(kT) 6 (t-kT) 


(la) 
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The Laplace transform of g*(t) is 


[g*(t)3 = G*(s) = I g(kT)e- kTs 
k=0 



(lb) 


( 2 ) 


Suppose that the analog transfer function G(s) is of the form 
m 

n 

i=l (s+ ai ) (3a) 

G(s)=K 

n 

n 

j=l (d+bj) 



1-10 


Where 

_3^ a complex zeroes of G(s) 

-b ■ complex poles of G(s) 
s » o + ju “ Laplace variable 
K = Constant 
n > m 


If there are no repeated poles in C(s) , then 


G(s) 


n 

E 

k»l 



Where is the residue at pole -b k * 


(3b) 


Since, 


Z[ 


s+u 


] = Z [ ae _ut ] - 


l-e- u Vl 


The standard z-transform of (3b) results in 


n 


G(z) = E Ri. 

k=l 


l-e~bkT z ~l 


(4) 


A third representation of interest is found by noting that multipli- 

n of g(t) and, E 6(t-kT) in the 
k-=o 

volution of G(s) and E e”^ s = 1 


cation of g(t) and, E 6(t-kT) in the time domain corresponds to con- 
k-=o 


k=o 

convolution the result is 


l-e -Ts 


in the frequency domain. After 


G*(s) - 1/2 g(CH-) + 1/T E G(s+jko) ) 

k=-« 3 


( 5 ) 
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where to =2ir/T. Hence it is apparent that G*(s) is periodic in to, the 
s 

sampling frequency. It is required that G(jto) = 0 for | to | > w s/2 as 
shown in Fig. 5 so that the frequency content of G(jto) will be preserved 

in the primary strip of G*(joo). If this relationship is preserved the 
envelope of G(joo) can be recovered from G*(ju). This phenomenon is 
known as frequency aliasing. 

The standard z-transform discussed above is best known of the 
sampled-data transforms. The theorems and tables of z-transforms can be 
found in any standard sampled-data text [1]. 

Hold Devices 


In the previous section we have seen the sampling process used to 
determine values of an input signal at discrete time intervals. The 
inverse process of data reconstruction from these sampled signals is 
accomplished by hold devices . Consider the problem of reconstructing the 
signal g(t) given samples spaced T seconds apart. If we expand g(t) 
in a Taylor's series 


g(t) = g(nT) + g' (nT) (t-nT) + ^ [ nT ) (t-nT) 2 + . . 


( 6 ) 


for nT <_ t <t-nT ,, 

where g'(t) = dg(t) 
dt 


The derivatives may be approximated by 


g' (nT) = £ (g(nT) - g(nT-T)) 


g"(nT) (g'(nT) -g'(nT-T)) 


( 7 ) 


» etc. 
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If equation (6) is truncated to just one term, this reconstruction is 
called a zero-order hold; if the first derivative is included, a first 

order hold; etc. 

Zero-Order Hold [1] 

The zero-order hold device in the mathematical model of Fig. 3 
accepts an impulse modulated input e Q *(t) and produces an output e Q (t) 
as shown in Fig. 2. The input e 0 *(t) may be expressed as 


e c *(t) - I e Q (kT) 6 <*-k.T) . 
k«o 

Its Laplace transform is 

E *(s) * E (z) = Z e (kT)z“ k 
o o •• o 

k=o 

The output (see Fig. 2) may be written 

OO 

e Q (t) = Z e 0 (kT) [u(t-kT)-u(t-kT-T)] . 
k=o 

Its Laplace transform is 

E q (s) - Z * o (kT)[ e ~ kTs ~ e~ kTs e' Ts j* 
k=o s s 

= Z e 0 (kT)e‘ kTs ( l-e-Ts . 
k=o s 

“ E 0 (z) j 


( 8 ) 
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The transfer function of the zero-order hold device is 


o ho (s) - VfL . « 


-Ts 


( 9 ) 


E o*< 8 > 


The frequency domain characteristics of G^ o (s) are shown in Fig. 6. 

Suppose that the sampling interval T is chosen very small. Then, 

e _T8 = 1-Ts +0(T 2 ) 

= 1-Ts 

Then, 

G ho (s) » 1 „ , ^ 1 ~ T g2. = T> (10) 

s 

This result is verified by noting in Fig. 6 that, for o><<u: 12 (or 

s 

T small), the magnitude function approaches T; also, the zero-order 
hold introduces phase lag which is linear with frequency into the 
system. 


Discrete Trans fer Functions [ 1 ] 


In the mathematical model of Fig. 1, the transfer function of the 
computing device is shown as 

E c (z) E *(s) 

« — - G(z). 

E^z) E^s) 

From Fig. 3 


E Q (s) = Ei*(s)G(s) . 
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Fig. 6. Gain and phase characteristics 
of a zero-order hold. 
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Hence, 

e *( t ) * l l e^ (kT)g(t-kT) ] l 6(t-jT) . 

k=c J = ° 

In this expression g(t-kT) may be replaced with g(jT-kT) due to the 
Dirac delta function. Next, let Jl-j-k and replace the summation 

index j with l as follows: 

00 00 

e *(t) “ E Z e. (kT) g(ZT) ,'5(t-ZT-kT) . 

° £=~V. 1 

Since g(y) 0 0 for y<0, the -k may be replaced by zero. The Laplace 
transform of e Q *(t) is 

E *(s) « [ Z eAkT)c KIS ][ Z gUT)e ] 

° k=o 4-0 

- E t *(8) G*(s). 

Hence, the descrete-time transfer function of the computing device in 
Fig. 3 is indeed 


B 0 *(s) 

.V<8> 


G*(s) 


G(z) . 


(ID 


Difference Equations [4] 

The discrete transfer function of equation (11) is, in general, 
the ratio of two polynomials in z"^ 


a 0 +a l z ~ 1+ ‘ * * +a n z_n E (z) 

G(z) = — 3 " — ° (12) 

l+b 1 z +..+b n z n E^z) 


where the coefficients a^ and bj are real numbers (can be zero). An 
equivalent expression for (12) is the equation 



1-17 


E 0 (z) = a 0 Ei(z) + a 1 z' 1 E i (z) + ... + a n z _ n Ei ( z ) 

-b 1 z“ 1 E 0 (z)- . .-b n z" n E 0 (z) . 

The infinite series for the z-transforms E (z) and E. (z) is now 

o 1 

substituted into the above equation and the coefficients of like 
powers of z _1 are equated, yielding 

e (kT) = a, e (kT)+a e, (kT-T)+. . .+a e, (kT-nT)-b, e (kT-T)-... 
o oi li ni lo 

-b e (kT-nT). (13) 

n o 

Note that z~l - e”Ts which represent a time delay of T second. Equation 
(13) may be programmed in the confuting device of Fig. 1* the variable 
ei(kT) is furnished by the A/D converter; delayed values of the input 
e^ (kT-nT) and delayed values of the output variable e 0 (kT-nT) are 
stored in the computing device. 

Equation (13) defines a programming scheme known as the direct 
form. A block diagram of this form appears in Fig. 7. 

Another programming scheme known as the canonical form is determined 
below. The transfer function is expressed as 

E (z) M(z) 

G(z) - 

M(z) E 1 (z) 
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where 

a Q + a 1 z”^+. . .+a tl z -n 
1 

l+biz - l+. . .+b n z" n (14) 

The time-domain equivalent expressions for (14) are 

m(kT )= e ± (kT) -b^kT-T) -. . . -b n m(kT-nT) 
e c (kT) = a 0 m(kT)+a 1 m(kT-T)+. . .+a n m(kT-nT) 


(15) 

( 16 ) 



Equations (15) and (16) are the difference equations to be used in 
the canonical programming form. A generalized block diagram of this 
form appears in Fig. 8. 



The standard z-transformation of an analog function in the s -plane 
may be considered to be a mapping from the s— plane to the z— plane under 
the rule 

z = e^ s . (17) 

See Fig. 9. The mapping illustrates that the region of stability in the 
s-plane (the left half-plane) corresponds to the interior of the unit 
circle in the z-plane. In fact, the primary strip in the s-plane maps 
onto the unit circle. All other strips also map to the unit circle 
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further illustrating the frequency aliasing problem of Fig. 5. Thus, 
transfer functions which are stable in the s-plane will also be stable 

after taking their standard z-transformation. 

Frequency Response 

The frequence response in the s-plane is evaluated by 

l G < s >l S =j0) 

/ G(s) s=ju> 0 to < Wg 

This corresponds to evaluating G(s) along the contour in the s-plane of 
Fig. 9a. Some upper cutoff frequence tog is shown for illustration. In 
the z-domain the contour follows the unit circle so that 

I D < z > I z = el“T 

/ D(z) z= ej“ T 0 < to < to B> 

is used to calculate the frequency response of a discrete transfer function. 
From equation (5) we see that D(z) is periodic in o) g so that in practice 

= u) s /2 and the contour traverses the top half of the unit circle. Hence 


db = 20 log | DCeJ^ 1 ) | 

<j> = /D(ej^ T ) 0 < m < io s /2 


( 18 ) 


will be used to calculate the frequency response of a digital filter. 


III. SAMPLED DATA TRANSFORMATIONS 


Sampled-data transformations are the techniques one uses to obtain 
numerical solutions to integral and differential equations. Any linear 
system's transfer function may be written as 


6 ( 8 ) 


= XM 

x(s) 


Y(s) = Laplace transform of the output 


X(s) = Laplace transform of the input. 


Alternately the relationship between input and output may be described 
as a differential or integral equation. Numerical methods may be 
employed to solve these equations; these methods approximate the integral 
and differential equations by difference equations. As we have seen pre- 
viously the difference equations may be represented by a discrete transfer 
function. The complete process is illustrated in Figure 10. 

Numerical Approximations 

Several numerical approximation techniques will now be presented, 
some for differentiation and some for integration. 

Backward Difference 

The backward difference is a simple technique which replaces the 
derivative of a function by 

1-23 
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_d_ (t) • y.Ct) - y(t - T) . 

dt T 

See Figure 11. 

In the Laplace domain 

sY(s) = Y.(s) - e~ sT Y(s) 

T 



T 


Hence , 


D(z) = G(s) 



T 


(19) 


Example . Find a discrete approximation for 


G(s) = — S_ 
s + a 

Y(s) = G(s) X(s) 


sY(s) + aY(s) = sX(s) 
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—[£- y(t) + ay(t) 


Now let 


y(t) = 


-k 


-- *(t - T) 
T 


Therefore 


Lt) - y(t - T) 
T 


+ ay(t) = ZL*) - ^ ~ T ) 


Evaluating at t = nT yields 


y(nT) = (x(nT) - x(nT - T) + y(nT - T)) 

1 + Ta 


Employing equations (12) and (13), 


D(z) = 


1 - z 


1 + Ta 


1 + Ta 
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An alternate solution employs equation (19) as follows 


D(z) = 


s + a 



aT + 1 - z -1 


D(z) = 

1 + aT 



1 + aT 


Forward Difference 


A similar numerical technique approximates 




4 y(t + T) - v f t) 


See Figure 11. 
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This represents the equivalent Laplace domain approximation 



T 


. z - 1 

S = 

T 
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D(z) 



D(z) = 

1 + (aT - 1) z” 1 . 

Rectangular Rule 

Suppose now we try some numerical approximations to integrals and 
compare results. 

Left Side Rule . Let us determine the numerical approximation for 

y(t) = } l x(t)dt . 
o 

Assume that the upper limit of the integral is t = nT. Hence 


y(nT) = / nT x(t)dt . 
o 


( 21 ) 


Figure 12a illustrates the rectangular rule using the left side of the 
rectangles. Hence 


n-1 

y(nT) = T E 


i=o 


x(iT) 


y(nT+T) = T E 
i=o 


n-1 

x(iT) = T E x(iT) + Tx(nT) 
i=o 


= y(nT) + Tx(nT) 
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Therefore using equations (12) and (13) 

Tz" 1 

D(z) = 

1-z -1 

= _T 

z-1 

Hence we have approximated the Integration transfer function 



which gives the same results as equation (20) for the forward difference. 

Right Side Rule . Figure 12b illustrates the use of the right side 
of the rectangle in approximating equation(21) . Therefore 

y(nT) = T E x(iT) 
i=l 

n+1 n 

y(nT + T) = T E x(iT) - T E x(iT) + Tx(nT + T) 
i=l i=l 

= y(nT) + T x(nT + T) 

Letting n = n - 1 

y(nT) = y(nT - T) + T x(nT) 

Employing equations (12) and (13) one finds 


D(z) = 


T 

1 - z“l 
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Hence, we have approximated the integrator 


I * _T 

s 1-z-l 


which yields the identical result of equation (19) for the backward 
difference. 


Trapezoidal Rule . 

The trapezoidal rule takes the average of the left and right side 
of the rectangles in Figure 12. Hence 


y (nT) - ~ tl E 1 [x(iT) + x(iT + T) ] 

L i=o 


X ii“X ti 

7 [ T E x(iT) + T E x(iT) ] 

2 i=o i=l 


Using the results of the rectangular rule, 


1 + z-- 1 

1-z" 1 


Thus we have approximated 


1^1 1 + z~ J 

s 2 1 - z”J 


This approximation is the familiar bilinear z - transform. 
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Simpson's Rule 

Simpson' 8 Rule evaluates equation (21) by the following formula 

y(nT) = 1 [x(0) + 4x(T) + 2 x(2T) + ... + 4x(nT - T) 

+ x(nT) ] 

But 

y(nT + 2T) = y(nT) + I [x(nT) + 4x(nT + T) + x(nT + 2T) ] 
Letting n = n - 2 and following equations (12) and (13) yields 

D(z) = T 1 + 4z -1 + z" 2 

3 1 - z=2 

Hence, we have approximated 

1^ * 1 + 4z - -*- + z“^ 

s 3 1 - z“^ 

Note that this formulation is valid only at even iterations (n even) . 
Impulse Invariance 

Suppose that we want to find a discrete equivalent filter for the 
Laplace transfer function G(s) . Further suppose that we desire the im- 
pulse response of the discrete equivalent to match that of the analog 
filter as shown in Figure 13. 

g(nT) = d(nT) . 

00 

D(z) = l d(nT)z~ n 
i=0 


Then 
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(a) Analog Filter 



(b) Digital Filter 


Figure 13. Impulse Invariance 
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= l g(nT)z~ n 
1=0 

= G(z) , 

which is the standard z-transform. Hence, for impulse invariance 

D(z) = Z[G(s) ] = G(z) 

the digital approximation is just the standard z-transform of G(s) . 


Impulse Invariant Integrator 

Let us find the digital equivalent of an analog integrator using 
impulse invariance and the models of Figure 14. We know that 


and that 



G hQ (s) = 1 - e~ Ts = T 


for small values of T. Hence 


V*> . 

X(z) 



and we have again approximated 



s 



Therefore, the backward difference, the right side rectangular rule, and 
the impulse invariant integrator all indicate equation (19) as their 
equivalent sampled-data transformation. 



X 


X(s) 


( 



( 


Figure 14. 








1-38 


Mapping Functions Summary 

As a result of our analysis of some elementary numerical approxima- 
tion techniques we have identified several sampled data mapping functions. 

Standard z-Transform 

The standard z-transform yields an impulse invariant filter the 
mapping function for this transformation is 

s = 1 In z . (22) 

T 

This mapping has been previously defined in Figure 9. 

Backward Difference 

The backward difference approximation for the solution of differ- 
ential equations provides the following mapping 

s = 1 - z \ . (23) 

T 

See Figure 15. Note that the region of stability in the s-plane maps 
into the right half plane z~ J > 1 of the z _1 plane. Since the region 
of instability in the z _1 plane is the interior of the unit circle, 
stable analog filters will always result in stable digital equivalents. 

In fact some unstable analog filters give stable digital ones. A major 
disadvantage of this mapping is seen in the frequency response contour. 
The ju> axis in the s-plane does not map to the unit circle in the z ^ 
plane (or the z-plane either). Hence, as we get farther from s = 0 or 
z = 1 the more degraded will be our desired frequency response. Thus, 
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we must decrease T (increase f 8 ) to improve this approximation. 

Forward Difference 

The forward difference approximation suggested the following 
mapping 

s = z - 1 . (24) 

T 

This mapping function is shown in Figure 16. Note that the left-half 
plane in the s domain maps to the region to the left of z - 1 in the z- 
plane. But the interior of the unit circle represents the stability re- 
gion in the z-plane. Consequently, some stable analog filters will 
give unstable digital ones. Unstable analog filters will also be un- 
stable digital ones under this mapping. Yet a further disadvantage is 
the same frequency contour encountered in Figure 15. Hence, this is an 
undesirable mapping. 

Bilinear z-Transformation 

The trapezoidal integration approximation led to the sampled data 
mapping 

s = 2 1 - z" 1 . (25) 

T 1 + z _J - 

This mapping is illustrated in Figure 17. Notice here that the entire 
left-half s-plane maps to the interior of the unit circle in the z- 
plane. Hence, all stable analog filters will result in stable digital 
ones. Also, the jw axis in the s— plane maps to the unit circle in the 
z-plane. However, the entire jto axis maps onto the unit circle which 
causes a mismatching of frequencies. This is a direct result of the 




2 



(b) z-plane 

Figure 17. Bilinear z-transform. 


1-43 


characteristic that for a digital filter 


z = 1 ->■ m = 0 
z = jl to = u > s/4 


z = -1 •> to 


(0 


S/ 2 


as required by equation (18). For the bilinear z-transform the frequen- 
cies in the z-plane (to^) are related to frequencies in the s-plane (^ A ) 
by 


or 


e ja) D T - 1 2j sin to D T 
JWA ” e j^T + 1 = ~ 

2 cos (Op! 

2 

top = 2_ tan 1 <o A 
T 


(26) 


See Figure 18. Correction for this frequency scale warping may be accom- 
plished by redesigning (prewarping) the critical frequencies of the de- 
sired transfer function G(s) before applying the bilinear z-transform. 

This transformation maps circles and straight lines in the s-plane 
to circles in the z-plane. It works well for frequency characteristics 
which are piecewise linear. It also insures that no frequency aliasing 
can occur in the transfer function frequency characteristic because the 


+j<*> axis does map into the upper half of the unit circle. Hence, the 
bilinear z-transform is quite popular. 
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Matched z-Transforms [2] 

The standard z-transform of equations (3b) and (4) required a 
partial fraction expansion of G(s) in order to complete the mapping 


s + u 


i a - uT ,-l 

l - e z 


For the purpose of simplifying the calculations, the matched z-transform 
maps the poles and zeroes (— bj and -aj of equation (3a)) to the z-plane 


as follows : 


s + a -*■ 1 — z" 1 e“ a ^ 


(27) 


Hence the matched z-transform of equation (3a) is 


G(z) = G(s) 


= K 


m 

n 

i = 1 


s + a ± - 1 - z -1 e _aiT 
s + bj = 1 - z“^e ^j^ 

(1 - z'V*! 1 ) 


n 

n 

j = i 


z- : e- a J T ) 


(28) 


where K is adjusted to give the desired gain at d.c. (z = 1) . This 
transform matches the poles and zeroes in the s and z planes. Note that 
the poles of this transform are identical with those of the standard z- 
transform but that the zeroes are different. Because of this difference, 
the matched z-transform may be used on nonbandlimited inputs. If G(s) 
has no zeroes, it is sometimes necessary to multiply (1 + z - -*-) N , N an in- 
teger, times the expression (28). 
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Other Transforms 

In general any transformation which maps the stable region of the 
s-plane into the stable region of the z-plane may be used. It is help- 
ful for the jw axis in the s-plane to map to the z-plane' s unit circle. 
Another important property is that rational functions G(s) should be 
transformed into rational functions D(z) so that the proper difference 
equations may be determined for realization. 

Simpson's Rule 

The Simpson's Rule approximation suggested that the mapping 

s = 3 1 ~ z~ 2 (29) 

T -1 -2 
1 + 4z + z 

be used as a transformation. The analysis of this mapping is left as 
an exercise for the reader. Please note that a second-order function 
G(s) will transform to a fourth-order D(z) . This is undesirable from a 
digital hardware viewpoint. 

(w.v) -Transform [14] 

In some applications, the system transfer function G(s,z,z a ) may be 
a function of s, z = e^ s , and z°, where 0 < a < 1. If all initial con- 
ditions are zero and 

w = 2. 1 - z 

T 1 + z -i 

v(a) = 1 - a(l - z' 1 ) + a (a - 1) (1 - z 

2 
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then for a system described by 

Y (s) = G(s,z,z a ) X(s) 
its z-transform will be 


Y (z) = G(w,z,v(a)) 


X(z) - x(0) 

1 + z -1 


If x(0) = 0, then 


D(z) = G(s,z,z a ) 


s = w 
z a = v(a) 

This completes the definition of the (w,v) transform. 


Example . Scott [15] has shown that a desirable phase lock loop 
has the transfer function 


G(s) 


10 


,-0.5 


s + lOz 

Using the (w,v) transform to find a digital equivalent if x(0) = 0 


D(z) = 


10 


s + lOz 


-0.5 


0 ' -1 
s = w = 2_ 1 - z 

T 1 + z -i 


z -0 - 5 = v(0.5) 


v(0.5) = 1 - 0.5(1 - z" 1 ) + 0.5 (-0.5) (1 - z 1 ) 2 


0.375 + . 75z -1 - . 125z -2 


D(z) = 


5T(1 +z -1 ) 


(1 + 1.875T) ' - (1 - 5.625T)z x + (3.125T)z - (.625T)z 


-2 



IV. DISCRETE STATE VARIABLES [5] 

An n fc k order discrete-time system is generally described by diff- 
erence equations. The difference equation description of the system 
dynamics may be alternately presented in vector matrix (state variable) 
form by the following set of first-order difference equations. 

x(kT + T) = Ax(kT) + Bu(kT) 

X(kT) = Cx(kT) + Du(kT) , (30) 

where x(kT) , u(kT) , and x(kT) are vectors of the discrete state variables, 
input variables, and output variables respectively. The symbol T is the 
sampling period and k is any non-negative integer. 

For the purpose of simplifying the notation, the sampling period 
T shall hereafter be omitted from the equation; thus, equation (30) be- 
comes 

jc(k + 1) = Ax(k) + Bu(k) 

£(k) = Cx(k) + Du(k) (31) 

The solution of equation (30) can be found by rewriting the first 
of equations (31) in standard z-transform notation: 

z X(z) = AX(z) + BU(z) + zx(0) , 

where x^(O) is a vector of the initial condition of the state variables. 
Solving for X(z) produces 

X(z) = (zl - A) _1 zx(0) + (zl - A) -1 BU(z) . (32) 
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The inverse z-transform (Z - *) of (zl - A) *z is 

Z -1 [(zI - A)” X z] = <Kk) = A k 

Therefore, the inverse z-transform of (32) yields 

, k-1 

x(k) = A K x( 0) + J A k -n Bu(n) . 
n=0 

This solution demonstrates that the present state of the system 
dependent upon the initial state x(0) and the system inputs u(n) 
the initial time (t = 0) to the present time (t = kT) . 


(33) 

x(k) is 
from 



V. CONVOLUTION 


In this section a review of continuous and discrete linear systems 
is presented. The equations are discussed in rapid succession. 

Continuous Linear Systems 

Figure 19a illustrates the conventional continuous system under 
consideration. Any linear system obeys superposition and is characterized 
by the impulse response g(t, £)* the response at time t due to an impulse 
at £. Hence 


y (t) = / x(£)g(t, OdK, 

—OO 

which is called the superposition integral, 
shift invariant then 

g(t, O = g(t - £) 

and equation (34) becomes 

OO 

y (t) = / x(£)g(t - OdC, (35) 

— OO 

t x(t)*g(t) 

which is the convolution integral. 

But, by definition of the Laplace transform 

00 

Y(s) = / y(t)e -st; dt . 

0 


(34) 

If the linear system is 
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Taking the Laplace transform of equation (35) produces 


Y(s) = G(s)X(s) , 

where 

00 

G(s) = / g(t)e~ St dt 
0 

is called the system transfer function; 

and G(jw) , the frequence response of the system. 

Discrete Linear Systems 

Figure 19b illustrates the conventional discrete-time system under 
consideration. The linear discrete system also obeys superposition and 
is also characterized by the impulse response d(nT, kT) , the response 
at time nT due to and discrete impulse 6 (kT) , where 


6 (kT) =1 if k = 0 

= 0 if k 4 0 


(36) 


and 


x(nT) = l x(kT)6 (nT - kT) 
k=0 


By superposition 


y (nT) = l x(kT)d(nT, kT) , (37) 

k=0 


which is called the superposition sum. If the system is shift invariant. 


then 
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d(nT, kT) = d(nT - kT) 
and equation (37) becomes 

oo 

y(nT) = l x(kT)d(nT - kT) , (38) 

k=0 

A x(nT)*d(nT) 

the convolution sum. But by definition of the standard z-transform 

OO 

Y (z) * l y(kT)z"\ 
k=0 

Taking the z-transform of equation (38) produces 
Y(z) = D(z)X(z) 

which is the identical result in equation (11). Hence, D(z) is called 
the discrete system transfer function and D(e^ a)T ) is called the frequency 
response (see equation (18)). 



VI. DISCRETE FOURIER TRANSOFRM 


This section examines the properties of continuous Fourier trans- 
forms and derives the discrete approximation. 

Continuous Fourier Transform 
The Fourier transformation may be defined by 

G(f ) - f g(t)e-J 2 ’ £t dt (39) 

—00 

and the inverse Fourier transform as 

g(t) = / G(f)e j2irft df . (40) 

—00 

The similarity between equations (39) and (40) is illustrated by the 
summary of Fourier transform pairs listed in Table 1. 

Another useful property of Fourier transforms is shown below: 

/“ |g(t)| 2 dt - /” |G(f)| 2 (41) 

_00 — CO 

This is known as the Fourier Integral Energy Theorem. 

Discrete Fourier Transform [6] 

Sampling Process 

In Figure 4, an impulse sampler was presented which sampled a 
signal for t 0. However, if the signal is zero for negative t, the 
following sampling function produces the same effect: 
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A(t) = l 6 (t - kT). (42) 

k=-c° 


From Table 1 the Fourier transform of A(t) is 

00 

FlA(t)] = ± l 6(f 
T k=-°° 

In order to verify this result we may note that F[A(t)] is periodic and 
may be expressed in a Fourier Series as 

£ -j 2irnTf 

F[A(t> ] = I c n e 
n=-“ 


where 


c 


n 


_ 1 _ 

T f 2 * F[A(t)]e j2lTnTf df 
2T 


c 

n 


_ 1 _ 

2T 

T / ^ df = 1 . 

± 1 

2T 


Hence 


F[A(t) ] = l 

n=-» 


-j 2irnTf 
e 


as expected from equation (42). 

Suppose the sampling function in equation (42) is multiplied by 
an input signal g(t) to produce the sampled signal 
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g*(t) = g(t)A(t) 

oo 

= l g(kT)<5 (t - kT). 

k=-00 

As illustrated in Figure 5, if B/2 tt, the highest frequency in the 
sampled signal is less than f g /2, then recovery of the original signal 
is possible with an ideal low-pass filter whose cutoff is f s . To 
recover the signal, g r (t) we multiply G (f) by a square window function 

V f) - 


G r (£) - G*(f)F fp (f) (43) 

Since multiplication of Fourier transforms represents convolution in the 
time domain 


g r (t) = g*(t)* 
r wt/T 


I g(kT)6 (t - kT) 
k=-“ 


, v - sitt ^(t - kT) 

g_(t) = l g(kT) 1 . 

^ i t r 


* 

sin(7rt/T) 


7Tt/T 


k=-00 


- kT) 


( 44 ) 


If f g > B/ir, the low pass filter is ideal, and the samples g(kT) are 
exact, then 


g r (t) = g(t) . 
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However, the samples are never exact, no signal is ever bandlimited, and 
no low-pass filter is ideal. Therefore, we can't exactly recover a 
sampled signal. 

DFT Derivation 

Now discrete versions for (39) and (40) will be determined. Define 
the following conditions for (39) : 

f g = sampling frequency 

T = l/f_ = sampling interval 

g(t) = 0 outside the interval t = [0,NT] 

N = an integer, the number of sample points 
G(f) is bandlimited to + f g /2. 

Please note that these conditions can never be completely satisfied. 

With the time-limited function g(t), G(f) cannot be bandlimited. In 
practice it can get quite small as |f| increases. However, since a 
function is never time and bandlimited, the time and frequency samples 
are corrupted by aliasing. 

Using the above conditions in equation (39) one obtains 
NT 

G(f) = / g(t)e" j2lTft dt . 

0 

Using the rectangular rule for numerical integration 
N-i rirn 

G(f ) = T l g(kT) e J 2irfkT . (45) 

k=0 
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The discrete version of equation (39) will compute samples of G(f) 
every Af = f g /N = 1/NT Hertz. Substituting f = nAf into equation (45) 

G(nAf) = T l g(kT)e -j2imAfkT 
k=0 


or 


N-l 

G(n/NT) = T l g(kT)e' J 
k=0 


(2ir/N)nk 


In another form 


G(n/NT) = T l g(kT)W _nk 
k=0 


(46) 


where 

W - e3 2 ’ /N . 

Since g(t) = 0 outside the interval 0<t<NT, one can construct 

a periodic function h(t) from g(t), with period NT: 

00 

h(t) = l g(t-mNT), (47) 

m=-°° 

which may be written 

t 00 

h(t) = / g (x ) l 6 (t - mNT - x)dx 

0 

00 

= g(t)* l 6 (t - mNT) 

m=-°o 
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where the * denotes convolution. Since G(f) is bandlimited to f g /2, 
so is H(f). Therefore, 


oo 

H(f ) = G(f) [1/NT l S (f - m/NT)] 
m=“°° 

- I Ht -m/m), 

IQS-00 


and 


CO 

H(f ) = l H' (m/NT) 6 (f -m/NT) , 
m=— 00 


(48) 


where H(f) is the continuous Fourier transform of h(t). The weighting 
function H' (m/NT) in equation (48) is defined below: 


H'(m/NT) - (l/NT)G(m/NT) . 


(49) 


Equation (46) may be inserted in (49), 


H" (m/NT) = (1/N) 


N-l 

l g(kT) W ' 
k=0 


From (47) , 


(50) 


g(kT) = h(kT) : k=0, N-l. 

Therefore (50) becomes 

N-l 

H' (m/NT) = (1/N) l h(kT) W _mk 
k=0 


( 51 ) 
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IDFT Derivation 


The inverse discrete Fourier transform is found by considering 
equation (40) 


00 

h(t) = / H(f)e j27rft df. 

— 00 


Under the conditions of the previous section 


1/2T 

h(t) = / H(f ) 


-1/2T 


eJ 2irft df 


(52) 


since H(f) is bandlimited to f _ /2 = 1/2T. 

s 

Substituting equation (48) into (52) and evaluating at t = kT, 


h(kT) = / l H'(m/NT)e^ 2irfkT 6'(f - m/NT)df 

-1/2T m=-» 

Since the integrand is periodic (1/T) in f, 

1/T °° 

h(kT) = / l H' (m/NT) e-* (27r/N)mk 6(f - m/NT) df 

0 m=-oo 

The limits of integration truncate the sum to 

h (kT) = l IT (m/NT) (53) 
m=0 

The prime is dropped in equations (51) and (53) for convenience, and the 
resulting relations are 
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DFT 

N-l 

H(m/NT) = (1/N) l h(kT)W _mk 
k=0 


(54) 


IDFT 

N-i . 

h(kT) * l HOn/NT)#* 
m=0 


(55) 


These equations define the discrete Fourier transform pair. This 
transform may be thought of as a mapping of N points in the time domain 
to N points in the frequency domain. 

DFT Pairs 

Although equations (54) and (55) are discrete approximations of 
(39) and (40), we can show that they form exact transform pairs. 

From equation (54) , 

H = | [W] h (56) 


where H and li are vectors constructed of the N samples in the frequency 
and time domains, and 



[i w-^ 
From equation (55) 
h = [W*] H 


1 

W -2(N-1) 

W-(N-l) 2 


(57) 


( 58 ) 
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where [W*] is the complex conjugate of [W] . In order to prove that the 
transform pairs are exact one must show that 


[W]" 1 = I [W*J . 

N 

This proof follows : 

Let 

[P] = [W][W*] 

then a general element of [P] is 

P £m = [1 W^... w" (N_1 )^ -1 )] 


(59) 


1 

y (m-1) 


|y(N-l) (m-1) 

- l + w" (£ " m) + ... + W - ( N_1 


Then all diagonal elements of [P] 


p ££ =l+l+...+l»N 


and the off diagonal elements 




. i - „ 


1 - W 


Hence 


[P] = N[I] 


and the exact relationship is proved. 
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V 


A s umma ry of DFT pairs is listed in Table 2. Several other 
interesting relations are displayed below: 

N-l 

h(0) = I H(m/NT) (60) 

m=0 

N-l 

H(0) = (1/N) l h(kT) (61) 

k=0 


and 

i N-l . N-l , 

» l |h(kT) | 2 = l | H (m/NT ) \ . (62) 

k=0 m=0 

This last relation is known as Parseval's Theorem. 


Fast Fourier Transform 

Calculation Time 

The fast Fourier transform (FFT) is a high speed technique for 
calculating the DFT. If the number of samples N may be written 

N " r i r 2 ,,,r n » r i an integer 

then 

[W] = fW 1 ][W 2 ]*“[W n ] (63) 

where [V^] is an N x N matrix with only r^N non-zero elements. The 
calculation of 

H - £ [W] h 
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TABLE 2. DFT Pairs 


h(kT) 

H(m/NT) 

Ah x (kT) + Bh 2 (kT) 

AH(m/NT) + BH(m/NT) 

h(kT - nT) 

W -mn H(m/NT) 

h 1 (kT)h 2 (kT) 

(m/NT) *H 2 (m/NT) 

h*(kT) 

H*(-m/NT) 

h(-kT) 

H(-m/NT) 

6 (kT) 

1/N 

6(kT - nT) 

(l/N)W mn 

, N-l 

± l h, 0£ + k)h 2 (£) 
N £=o 

H 1 (m/NT)H 2 (-m/NT) 
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2 

requires N operations of complex multiplication, whereas 

H = I [W 1 ][W 2 ]---[W n ] h (64) 

^ 

r n N operations 

requires (r^ + r 2 + ••• + r n ) N operations. For the special case r^ = 2, 

N ■ 2 n , the total number of operations is 

#oper = (2 + 2 + *•• + 2)N 
= 2nN 

= 2Nlog 2 N . (65) 

Example . Compare the time to calculate the DFT and FFT of a sequence 
of 1024 samples of a time function given that a typical computer cal- 
culates a complex multiplication in about 40ys 

DFT: 

Calc. Time = N^(40ys) 

= 1.053 x 10 6 x 40 x 10 -6 
= 42.1 seconds 

FFT: 

Calc. Time = 2Nlog^N(40ys) 

= (2048) (10) (40 x 10 -6 ) 

= . 82 seconds 
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FFT Derivation 

Since the DFT is a linear operation and N = 2°, we may break equation 
(54) into two functions, the even samples and the odd ones: 


H(m/NT) 


(N/2)-l „ , _ m 

I I [h(2kT)W 2mk + h(2kT +T)V~ m \ ] 

N k=0 


1 

N 


(N/2)-l -m (N/2)-l 

l h(2kT)W -2 ^ + l h(2kT + T)W _2mk 

k=0 N k=0 


or 


H (m/NT) = DFT [h(2kT) ] + W -m DFT[h(2kT + T) ] 


( 66 ) 


for m = 0, y - 1. But 


— m — 


N 


N 


W 


2 = W m w" 2 = -W 


-m 


N 


w -2(m + 7) = w -2V N - W' 2 " . 


Therefore, the remaining samples may be determined by 

(N/2)-l _ 2 , 

H(m/NT + 1/2T) l h(2kT)W _ 

N k=0 


W 


-m (N/2)-l 


N 


l h (2kT + T)W 


-2mk 


k=0 


or 


H(m/NT + 1/2T) = DFT[h(2kT) ] - W -m DFT[h(2kT + T) ] 


(67) 
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for m = 0, y - 1. The above equations may be successively applied in 
order to achieve the maximum reduction in computation time indicated 
by equation (65). The technique of dividing the time samples into 
even and odd parts is sometimes called "decimation in time." The FFT 
of eight-points is illustrated in Figure 20. 

Note that the time samples are entered in "bit reverse" order: 




bit 



binary 

reversal 


0 

000 

000 

0 

1 

001 

100 

4 

2 

010 

010 

2 

3 

Oil 

110 

6 

4 

100 

001 

1 

5 

101 

101 

5 

6 

110 

Oil 

3 

7 

111 

111 

7 


I FFT Derivation 

Since the IDFT is a linear operation and N = 2 n , we may separate 
equation (55) into two functions, the even samples and the odd ones: 


(N/2)-l 

h(kT) = l [H(2m/NT)W + H(2m + 1/NT)W W k ] 
m=0 


(N/2)-l (N/2)-l . , 

l H(2m/NT)W ZmK + VT £ H(2m + l/NT)bT mk 
m=0 m=0 


or 


h(kT) = IDFT[H(2m/NT) ] + W k IDFT[H(2m + 1/NT)], (68) 


for k = 0, y - 1. But 
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Figure 20, FFT for eight samples. 
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N N 

/ + 2 , „k„7 . 

w 2(k + j) a w 2k„N „ w 2k _ 


Therefore the remaining samples may be calculated by 


, NT, 
h(kT + -y) = 


(N/2)-l . , , (N/2)-l 

l H(2m/NT)W^ mlc - W k £ H(2m + 1/NT)W 
m=0 - m=0 


,2mk 


or 

h(kT + - IDFT[H(2m/NT) ] - W k IDFT[H(2m + 1/NT)] (69) 

N 

for k = 0, — - 1. The repetition of this process yields the algorithm 
for the IFFT. The above derivation is sometimes called "decimation in 
frequency. " 

Equations (66) through (69) suggest an algorithm for calculating 
the IFFT as shown in Figure 21. In this figure the equations (68) 
and (69) are employed at each stage of the transformation of eight 
frequency samples into eight time samples. Note that the frequency 
samples are again inserted in "bit reverals" order. 

The reader will please note the similarity of the Figures 20 and 
21. The basic element is sometimes called a "butterfly" as shown in 
Figure 22. The gains on a few multipliers are different. This structure 
suggests the mechanication of FFT hardware. 
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VII. RANDOM PROCESSES 


In the analysis and synthesis of digital filters one frequently 
encounters signals which are random in nature that must be examined 
with special techniques. 


Continuous Processes [6] 

If G(f) is the Fourier transform of a continuous signal g(t), the 
power density function or power spectrum may be defined as 



The auto correlation function 




gg 


(t) = lim 
A-*» 


1 

A 


A 

/ g(t)g(t + t ) dt. 
0 


These two functions form a Fourier transform pair 

oo 

4'gg( f ) = / ’i'gg('r) e J2irfT dt 

Both the power spectrum and the auto correlation function are real and 
symmetric. 

Discrete Processes [16] 

For the discrete case, the cross-correlation function is first 
defined 
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1 N-l 

ip (kT) - lim - l y(nT)x(nT + kT) 
xy N-~» n-0 


When both functions are the same (x = y) , the cross-correlation becomes 
the auto -correlation 

i N-l 

^(kT) « Hm n l x(nT)x(nT + kT). 

N-h» n=0 

The auto-correlation function evaluated at k = 0 yields 

i N ” 1 ? 

il* (0) = lim - l x (nT) 

N-**> N n=0 

= x 2 (kT) 

the mean squared value of the signal x(kT). 

Since the power spectrum is the Fourier transform of the auto- 
correlation^we see from Table 2 that 

4^ (m/NT) * X (m/NT) X (-m/NT) 

if the signal x(kT) is time limited in the interval [0, NT]. Therefore 
4^ (m/NT) = | X (m/NT) | 2 . 

The Fourier transform of the discrete cross-correlation function is 


4* (m/NT) = X (-m/NT) Y (m/NT) 

Ay 


and is sometimes called cross power spectrum or cross-periodogram. 
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Lastly, consider the discrete system of Figure 13b with a random 
input whose power spectrum is known. We may find the mean squared 
value of the filter's output by the following: 

y 2(nT) = “y / 'f xx (z)D(z)D(l/z)dz/z 

where 

T = the unit circle 
00 

V (z) - l *xx( kT > z " k 

k=0 


= power spectrum 
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I. DIGITAL FILTER CATEGORIES 


The generalized transfer function of a digital filter has been 
shown to be the ratio of two polynomials in z. Generally the coefficients 
of the polynominals are real numbers which must be determined in some 
manner to force the digital filter transfer characteristics to meet 
some criteria. The manner in which the coefficients are found as well 
as other considerations (for example, implement ion details) allow us 
to categorize digital filters into several classifications. In this 
chapter we examine non-recursive, recursive, and several other major 
categories for digital filters. 

Non-Recursive Filters 

General 

Non-recursive digital filters are those whose transfer function 
can be written as 

D(z) = £ a i z“ i . ( 1 - 1 ) 

i=0 

Non- recursive filters have no feedback terms, and hence they have a 
finite impulse response. They are sometimes called transversal filters, 
a name used for delay line filters in radar moving-target-indicator 
applications. 
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Much emphasis has been placed on these filters in the literature 
and several design techniques will now be illustrated. 

Finite Impulse Response Filters [1] 

Finite impulse response (FIR) filters satisfy equation (1-1). 
Consider the first order filter 


H(z) = 


1 - az 


-1 


a <1 


= 1 + az“l + (az -1 )^ + ••• 


= l (az 1 ) /£ . 

1=0 


Suppose we truncate the series to M terms to produce the FIR below 


M-l , o 
H m (z) = l (az" 1 ) 

1=0 


( 1 - 2 ) 


Also , 


. , -1. M 

H m (z) = (az > 


- a M - 


1 - az 


-1 


z M-1 (z - a) 


(1-3) 


This is another way of expressing the FIR as filter with feedback. 
Fig. 1 illustrates the z-plane pole-zero locations for both H(z) 
and H m (z) for M = 8. Fig. 2a shows the implementation of (1-2) ; 
Fig. 2b, (1-3). 
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Fast Convolution [2] 

Fast convolution is a technique which employs the FFT and IFFT 
to determine the filter output response (see Fig. 3a). Direct convolu- 
tion is expressed by 

N-l 

y (kT) = l h(£T)x(kT - £l) . (1-4) 

1=0 

2 

To calculate N output points, this requires N real multiplications. 

For fast convolution 

y(kT) = IFFT{H(— ) *FFT(x(kT)) } . (1-5) 

NT 

Here the FFT and IFFT require 2Nlog2N operations each, while the multi- 
plication requires N operations. This totals 

// operations = N(4 log 2 N + 1). 

If each operation (a complex multiplication) is assumed to take approxi- 
mately 4 real multiplications, the result is 

// multiplications : 16 N log 2 N. (1-6) 

Suppose N = 1024, then ~ 10^ and 16 N log 2 N - 1.6 x 10^. Hence, 
for large numbers of output points, the fast convolution technique 


is faster than direct convolution. 
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The reader must be careful in using the fast convolution technique 
because the results can be misleading. Consider the convolution of 
the analog signals in Fig. 3b. If we sample these signals and use 
the FFT and IFFT, we are convolving the periodic functions shown in 
Fig. 3c. Hence the output y(kT) can differ greatly from the desired 
sequence. 

In order to improve the results one may add zeroes into the input 
and transfer function sample sequences as shown in Fig. 4. Note the 
improved output response. However, in adding zeroes we have increased 
the calculation time unless we modify the FFT algorithm. 

Linear Phase Filters [3] 

A linear phase filter is an FIR filter with exact linear phase. 

They may be used to approximate an arbitrary magnitude frequency 
response without causing phase errors. The linear phase filter is 
good for standard lowpass » bandpass, and highpass filters. 

If the number of sample points in a FIR filter is 

N = 2t + 1 , 

then linear phase with delay t is realized if and only if 


h(kT) = h(NT - T - kT) 


(1-7) 
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Hence, for N even 

1. There is no unique peak in h(kT) 

2. h(kT) = h(NT - T - kT) k = 0, (If) - 1 

3. The center of symmetry is between (Ip) and (Ip) - 1. 

4. The delay is 

t = N ~ 1 
2 

the center of symmetry. 

For N odd 

1. There is a unique peak in h(kT) at (N - l)/2. 

2. h(kT) = h(NT - T - kT) k = 0, (N - l)/2 

3. The center of symmetry is at (N - l)/2 

4. The delay is 

_ N - 1 
T 2 

the center of symmetry. 

If the above conditions are met, the frequency samples H(^.) will be 
given by 

h< st> ’ eJe ” 

where, for N even 
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and for N odd 

6 = - ^ mt m = 0, (N - l)/2 

N (1-9) 

e m = -Y < N ” m ) T m ■ (j) + 1. N . 

and 

1 f s 

H(— ) = H(— ) = 0 . (1-10) 

2T 2 

This concludes our brief description of linear phase filters. 

Frequency Sampling Filters [31 

The term frequency sampling filters refers to a class of digital 
filters specified by sample points in the frequency domain and imple- 
mented in the manner of Fig. 5. Many techniques have been suggested 
for choosing the sample points 

H( OT } = l H ml ej9m m = 0, N - 1 (1-11) 

including optimization techniques which adjust the points in the 
transition region to give a good ripple between sample points. For 
real impulse response filters 



«.*#■*■* 

Pass 

Band 



JL 


0 


x(kT) 



Fig. 5 


2 -] 



NT 



(b) 


Frequency Sampling Filter 



2-12 



6 


m 


0 


N-m 


By the IDFT 


( 1 - 12 ) 


N-l 

h(kT) = l H(— ) 
m=0 NT 


e j (2Tr/N)mk 


k = 0, N - 1 


and 


N-l 

H(z) = l h(kT)z“ k 
k=0 


where 


H (z) 


= <) 


z = e j27Tm/N 


Then 


H(z) 


N-l 


I 

k=0 


r N-i 


l 

m=0 


H(-2L) 

NT 




-k 

z 


or 


H(z) = (1 - z" N ) 


N-l 


l 

m=0 


H*) 


1 - z'^w® 


(1-13) 
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where 


W = 


e j2ir/N 


Equation (1-13) is the motivation behind the frequency sampling 
implementation illustrated in Fig. 5b. 

Windowing Filters [4] 

In equation (1-2) a truncation was performed (a fairly drastic 
measure) to produce a FIR filter from an infinite impulse response 
function. Windowing is the process of orderly termination of an 
infinite series by truncating the series and adjusting the remaining 
terms to mask the truncation effects. The transfer function for the 
FIR is given in equation (1-1); its output response is 

m 

y(kT) = l a.x(kT - iT). 
i=0 

Briefly stated, the problem is to find the coefficients a.^ of the 
FIR filter H(z) such that 


H(e j “ T ) .= F^" 1 ) 


where F(e^ b> '*’) is some specified desired frequency response. The 
design procedure is outlined below: 
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1. From F(e-^ w '*'), use the IFFT algorithm to find f (kT) . 

2. Multiply f(kT) by a window function w(kT) , or 

h(kT) = f (kT)w(kT) (1-15) 

The process is outlined in Fig. 6. Multiplication in the time domain 
is convolution in the frequency domain, and hence 

H(f ) = F(f ) *W(f ) . 


The window function shown is a rectangular one which duplicates the 
truncation process. Notice the ringing effect in Fig. 6c. The side- 
lobes for this window are about 20%. Fig. 7 illustrates two other 
windows. The triangular one reduced the sidelobes to about 4%. The 
raised cosine window is the best one shown. Its function is 


w(t) = a + (1 - a) cos 


TTt 

T 


(1-16) 


If a = 0.50 it is called a Hamming window. The optimal value of a is 
about 0.54. This value yields the Hanning window and reduces the 
sidelobes to about 1%. 


Moving Average Filter [5] 


A moving average filter is a FIR filter which calculates the 


average of the N most recent observations of the input: 



(a) Desired Frequency Response 


W(w) 




(c) Windowing Filter 


Fig. 6. Windowing Filter Construction 
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w(t) 



Fig. 7. Window Functions 
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y(kT) - i l x(kT - ll) 

N l- 0 


and 


, N-l 

H(z) = £ l z 
1=0 


-l 


In another form 


H(z) 



(1-17) 


Least Mean-Square Digital Filters [6,7] 

Assume that the filter Input Is x(nT), a random signal whose 
autocorrelation R X x( t ) I s known, and that the crosscorrelation 
Rd x (t) of this input and a desired output d(nT) is also specified. 
Let the impulse response of the filter be g(kT) , and its output, 
z(kT). Allowing a shifted time scale, 

N 

z(nT) = l x(nT - KT)g(kT) . 
k=-M 

Define the signal D to be expected value of the difference in the 
actual and desired filter outputs squared: 


D = E[d(nT) - z(nT)] 2 . 
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By definition 


R (t) = E {x ( t + t)x(t)} 

R dx ( c ) = E t d (t + T)x(T)} (1-18) 

R dd (t) = E {d (t + T)d(x)} 


Substitution of these relations into D yields 


N 

D " R dd (0) ' 2 I R dx (nT)y(nT) 
n=-M 

(1-19) 

N N 

+ l l R* x (kT - nT)g(kT)g(nT). 

k— M n— M 


The purpose of the least-mean squares filter is to minimize D by 
choosing g(nT). Hence, if one takes the partial derivative of D with 
respect to g(nT) and sets the result to zero, the following solution 
is generated 


N 

R dx (nT) = l Rx X (kT - nT)g(kT) . -M <_ n < N (1-20) 

k— M 


In equation (1-20), all quantities are known except g(kT). Hence, 
the filter weights may be calculated from (1-20). Least mean-squares 
filters are sometimes called digital Wiener filters. 
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Least Squares Polynomial Moving Arc Filter [5] 

The problem here is to solve for the coefficients a^ of a polynomial 
to best fit the input data y(t^) in a least squares sense. Each input 
point is approximated by 


y(t i ) = aQ + a-j^ + ^2^ + + a d c i 


- I 


k*0 


Vi 


If the input samples are evenly spaced, t^ * iT and 


y(iT) = l a, (iT)' 
k=0 


For n input samples define 


n 

s- l 

i=0 


l a. (iT) k - y (iT) 
k=0 


In order to minimize S by choosing a^, one may take the partial 
derative 



l a. (iT)' 6 - y(iT) 
k=0 



= 0 


This expression reduces to 
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d d , d „ 
l l a k (iT) {lT) Z = l y(iT) (£l) 
k=0 1=0 1=0 


Written in matrix form 


CA = B 


or 


A ■ c 1b - (1-21) 

In (1-21) the matrix C _1 represents the filter itself (whose coefficients 
are precomputed) and B represents the system input. The output is A 
which represents the polynomial coefficients a k . 

Another form of polynomial filtering termed exponential filtering 
allows the polynomial to grow by one term as each new input occurs. 

Such schemes are called "growing memory" filters. 

Digital Inverse Filtering [8 .91 

Digital inverse filtering is a special case of least mean-square 
filtering as described in equation (1-20). Suppose that the desired 
filter output is 

d(kT) = 1, 0, 0, 0, 

the discrete impulse function. Hence the crosscorrelation 



2-21 


R dx (nT) = x(0) » 0, 0, 0, ••• 


which can be scaled to unity (x(0) = 1). Equation (1-20) with M = 
then becomes 


* r 0 r 2 ••• r N 


'so' 


r l ’ 

r l r 0 r l r N-l 

• • • 


»1 

• 

= 

0 

• 

r N r N-l r N-2 r 0 


• 

®N J 


• 

0 


where 


r i = r -i " *«<«> * 


If the filter coefficients are used in an FIR filter 


H(z) = l g^” 1 
i=0 


( 1 - 22 ) 


and the random signal x(nT) is applied to the input, the output will 
be a digital impulse function. Therefore, H(z) is said to be an 
inverse digital filter. The calculations involved are shown in Fig. 8. 
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*xx< iT > 


d(kT) 



g, 


(a) Filter Design 


x(kT) 




z(kT) = d(kT) 


(b) Filter Application 


Fig. 8. Digital Inverse Filtering 
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Recursive Digital Filters 

General 

A recursive digital filter is a filter with feedback which, in 
general, has an infinite impulse response. Its transfer function is 


n 



H(z) = (1-23) 

1 + l biz" 1 
i=l 

where at least one a^ and b^ is not zero. 

Recursive filters generally require fewer terms (lower order) 
than a non-recursive filter with similar characteristics. Higher 
order recursive filters are usually factored into second order stages 
which are either cascaded or paralleled. 

Block Recursion [10] 

One technique for implementing a desired recursive digital filter 
of the form 


Hp(z) = 


1 

D(z) 


(1-24) 


is called block recursion and is shown in Fig. 9. The implementation 
in Fig. 9b is 
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H (z) 

v z) = V - - 

1 - Z ^(z) 


(1-25) 


But the desired filter is 


V 2> • T 


n (l - z z" 1 ) 
i=l 


m a 

1 — 
i=l 1 - z^z 


-1 


where z^ are the poles of the function Hp(z). 

The finite impulse response filter H^/z) is found by truncating each 


component of Hp(z) to M terms, or 


H m (z) = 


® Sj^ll - (ZjZ -1 ) 11 ] 
i-1 1 - z ± «-l 


= H p (z) - z-M l 

V 1*1 1 - ZjZ 


-1 


1 - z -M Q(z) 
D(z) 


(1-26) 


where 


m a.zVD(z) 

Q(z> " l 11 


i-1 1 - ZjZ -1 


Thus, 


H m (z)D(z) = 1 - z _M Q(z) 
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and 


Q(z) = z M (l - H m (z)D(z)) 


( 1 - 27 ) 


where Q(z) is a polynomial of order M-l. From the above relations it 
is clear that if G(z) in Hg(z) is chosen as Q(z), then 


H b (z) 


" H p (z) 


G(z) = Q(z) 


( 1 - 28 ) 


and the block recursive implementation exactly produces H (z) , the 

P 

desired filter. Thus, we have shown that a recursive filter can be 
implemented using one FIR Hjj(z) in the feed forward path and one FIR 
G(z) in the feedback path, where 


H^(z) = truncated version of Hp(z) 

M < X - 29 > 

G(z) = z (1 - H M (z)D(z) ) . 

Some researchers have used the FFT to implement the two FIR filters 
[ 11 - 13 ]. 

Example . Consider the filter 

1 = _1 

1 + az - -*- + bz - ^ D(z) 


H p (z) - 



2-27 


and let M = 3 


1 2 . 1 ~ az + (a 2 ~ b)z~ 2 

1 + az -i + bz /~1 

1 + az~^ + bz - ^ 

- az~^ - bz -2 

- az~^ - a 2 z~ 2 - abz~ 3 

(a 2 -b)z“^ + abz -3 


Hence 


H 3 (z) = 1 - az“^ + (a 2 -b)z- 2 
G(z) = z 3 (1 - H 3 (z)D(z>) 

= (2ab - a 3 ) + (b^ - a^b)z”^. 


One can check the impulse response of Hg(z) by dividing the 
denominator into the numerator and comparing it with Hp(z). 

Flat Group Delay Digital Filters [141 

In order to achieve a linear-phase digital filter one must choose 
a non- recursive structure. However, when the order of the non- recursive 
filter is unacceptably larger, one is led to approximate the linear 
phase using a recursive filter design whose error norm is the maximally 
flat criteria. 

Consider the recursive filter 
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H (z) = 


1 + l a i 

i=l 


0 


1 + 7 a.z 

i=l 1 


-i 


whose d.c. gain is unity. The phase response (T = 1) is given by 


tan 


-1 


n 


l ai sin 1 

i=0 


1 0 ) 


n 


£ a^cos i a) 
i=0 


= $(oj) 


(] 


The ideal phase (—to t) , where x is the desired delay is approximated 
by minimizing 


6 (to) = -IDT - $(<!)) 


0 


or 


e(oj) = tan(<5(io)) = - tan wt - tan ($(w)) . 

The procedure is to make e(oo) vanish at d.c., together with its 
derivatives up to some order depending on n. 

The solution yields the filter 


-30) 


-31) 


-32) 
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H(z) 


2n! 1 

n! 2n 

n (2 t + i) 

i=n-KL 


1 


n i 

l (-D k 

k»0 


n 



2t + i 
2 t + k + i 



which is stable for all finite positive values of t . 


( 1 - 33 ) 
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Advanced Topics 

In addition to the simple division of digital filters into recur- 
sive or non-recursive categories, there are many other ways of identifying 
their characteristics. 


Complex Digital Filters [15] 

A complex digital filter has a complex input x(nT) , a complex 
output y(nT), and a complex transfer function H(z). An example is 
shown in Fig. 10. A lowpass envelope is centered at f by replacing 
z by e'j^c^ z = yz in H(z). 


H(z) = 


l a Z z 

1=0 


-i 


i + y bpz 

i=i 1 


-i 


Hence 


Hg (z) 



(1-34) 


Complex digital filters have application in communication and 
information theory, signal detection, randomly time-invariant channels, 
etc. In one application they are used to generate the Hilbert transform 
of a real signal x(nT). 
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,H s (f) 



Fig. 10. Complex Bandpass Filter 
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Randomly Sampled Filters [16] 

A randomly sampled digital filter H(z) takes input samples of the 
analog input signal x(t) at some random time 

nT 1 t n i nT + T 

and stores them in an input buffer. The numbers x(t n ) are then fed 
to the filter hardware as evenly spaced samples x(nT). Hence, the 
direct convolution of x(nT) and h(nT) produces the output y(nT) 
which is interpreted as y(t n ). The question arises what errors are 
generated by the random sampling? 

Define 


t„ - (n + Z n )T 

- A < Z < A 

— n — 



0 < a < 1 

where is a random variable. Also define 
n 

x(t n ) = x n = x(n + V 


(1-35) 


x(nT) = x . 

n 


(1-36) 
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a 

If we expand x n in a Taylor series about 2^ 

^ * *n + Vn + 1/2i n z n + "' 

• d 

where x R ■■ Xp, and define an input error C n 
?n = *n ’ x n - %Z n + 1/2 \ Z n + 

The output error due to random sampling is defined as 

% */• - »» ‘ j 0 h »-i 5 i 

for a non-recursive filter. Other pertinent relations are 

v 2 = ECZ 2 ) 
n 4 = E(Z 4 ) 

E(C n ) - 1/2^2 + ••• 

E(C 2 ) = x 2 v 2 + (l/4x 2 + l/3x^x n )n 4 + “ * 

E(e n ) * I h n _ E(q) 

- 1 — - - 

E(e 2 ) = E 2 (e n ) + [ h^varCS.) . 

i-0 
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The frequency response error is determined for sinusoid inputs 


x (t) = cos iot 0 <_ to £ Nf 


where N f is the Nyquist frequency ir/T. The expected values of the 
output steady state errors are 

E( e n ) ss = (-l/2io^v2)H( e J a) ) cos nio 


_ / 2 \ _2, N , / ? .2 \ /m 2m2 ,A,4 ,A,4 

E < e „ >». ■ E < e „> 8S + Z *>i 

i=0 


<o*v* <o H v H oTv 




2 24 

7io^v^ lo^v^ 

~ 24 8~ 


cos2nio 


The physical interpretation of the results is shown in Fig. 11, where 


EO^e^)) = (1 - l/2io 2 v 2 )H(e^) 


(1-38) 


In random sampling only the expected amplitude response is distorted 
while the expected phase response is unchanged. The noise to signal 
ratio for noise generated by random sampling is approximated by 


NSR = 10 log 1Q 


£ h? u)2v 2 
i=0 


|H(e^) 


iuK I 2 


(1-39) 
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(a) Exact Model 



(b) Approximate Model 


Pig. 11. Randomly Sampled Filter 
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Random sampling finds wide application in time-sharing filters, 
radar filters, and faulty samplers - all samples are faulty to some 
extent . 

Example . [9] 


1 - . 9z -1 

and Z n has a rectangular distribution with a = 0.1, or 10% jitter 
in the input sampler. The curves of Fig. 12 illustrate that as 
frequency increases, the noise component increases making the 
filter unusable above w/N^ =0.3 

Multirate Digital Filtering [171 

A multirate digital filter is one in which the samplers for the 
input and output are operating at different rates, one usually being 
an integral multiple of the other. Much analysis of multirate sampled 
data control systems has been treated in the open literature. Here 
we examine three configurations of multirate filters demonstrated 
in Fig. 13. Solutions for the sampled output frequency responses are 

00 

W**(ju>) = _L l G(jm + j |m) R*(joj + j Ip.) 

KT n=-« Ki 
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R(s) 


/ R*(s) 
T 

fast 



W(s) / W*(s) / W**(s) 

T KT 

slow 


(a) Typical Multirate Filter 


R(s) 



(b) Digital Prefilter 


R(s) 



slow 


(c) Analog Prefilter 


Y**(s) 


Fig. 13. Multirate Digital Filters 
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00 - 

X**(juO .1] G**(jw)H(ju) + ^ ) 
KT n=-®L 


Y**(jw) 


00 

4 I < 

KT n— « . 


G**(ju))F(jui + ) 

KT 


R*0» + -iff ) 


R(Ju + ) 

KT 


(1-40) 


These expressions simplify greatly if the filter function G(s) is 
band limited 


I G(joj) | = 0 


M > — 

KT 


(1-41) 


The functions H(s) and G(s) represent prefilters used to band limit 
the input signal r(t) to prevent frequency aliasing. 


Two Dimensional Digital Filters [18,19] 

Two dimensional digital filters are used in digital image 
processing. They are used to transform characteristics in photographs 
or CRT images. The transformation is described by 


H(z^ ) Z£) 


l l 

i=0 n=l 

! I 

1=0 Ti = 1 


m n 
a mn z l z 2 


i _m_n 
b mn z l z 2 


A(z l> z 2 ) 
B(z^ > ^2) 


(1-42) 


and 
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-Si A 

Zi ** e ■*- 


-s B 

z 2 = e 2 


where and are Laplace variables; A and B are the sampling intervals 

in the x and y planar coordinates of the image being processed. The 
two dimensional filter may also be expressed as 

00 00 

H(zi»z 2 > “II h mn z l z 2 (1-43) 

m=0 n=0 

where h mn is the impulse response. 

Let us consider the stability of a two dimensional digital filter. 
For a stable filter 


I I I hnm I < 00 • (1-44) 

m=0 n=0 

A two dimensional digital filter is stable if and only if no value of 
zj and Z 2 exist such that 

B(z^,Z 2 ) = 0, and 
I Z 1 1 — 1 » and 

I z 2 I 1 !• 
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Equivalent conditions are listed below: H(z 1 ,z 2 ) is stable if and 

only if 

1) The map B(zj,z 2 ) * 0 of the unit circle |zjJ = 0 to the z 2 

plane is outside the unit circle |z 2 | ■ 1. 

2) No point in | z^ | <_ 1 maps into z ^ - 0; or z 2 - 0 maps 

outside the unit circle in the z^ plane. 

Example . Given the two-dimensional filter 

H(z 1 ,z 2 ) = - 

1 + azj + bz2 

we may set the denominator to zero 
B(z^,z 2 ) = 1 + az^ + bz 2 = 0 
to determine the following map 



Condition 1 is shown in Fig. 14. 
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Fig. 14. Stability in Two-Dimensional 
Digital Filters. 
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or 


a| + | b | < 1 


for stability. 

i 

Condition 2 checks the point z 2 = 0 in the plane: 



but Zj must lie outside the unit circle, so 



a| < 1 


which is included in condition 1. Therefore, the example filter 
is stable if the sum of the magnitudes of the coefficients is 
less than one. 

Non-recursive two-dimensional digital filters may also be designed 
using windows, just as their one dimensional brothers. If w^(x) is a 
good one-dimensional window, then 
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w 2 (x,y) = w x ( / x 2 + y 2 ) 

will be a good two-dimensional window function. 

Example . Consider the one-dimensional window 

w(x) » 1 - |x| |x| < 1 

= 0 |x| > 1 

Then 

W(u) = sin 2 (to/2 )/ (io 2 /4) 


which has sidelobes of about 4%. 
The two-dimensional counterpart is 


w 2 (x,y) = 1 - »^ x 2 + y2 


x 2 + y 2 | <_ 1 


= 0 


x 2 + y 2 | > 1 


Then, in the frequency domain 


W 2 (o) 1 ,oj 2 ) = 2ir [p“ 3 / J 0 (t)dt-p 2 J o (p)] 

0 




to 2 + to? 


(1-45) 


which has sidelobes of only 2%. 
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The reader is referred to the open literature where much two- 
dimensional digital filtering theory is reported. 

Adaptive Digital Filters [20] 

A major advantage which digital filters hold over analog ones is 
the ease in which a digital filter's coefficients may be changed while 
the filter is processing data. Adaptive digital filters change their 
coefficients to minimize some specified criteria. An example non-recursive 
adaptive digital filter is depicted in Fig. 15a. The filter output is 

K ( 1 ) 

d(iT) = l g} x(iT - kT) (1-46) 

k=0 k 

where are time varying coefficients calculated as shown in Fig. 15b. 

g£ L+1) = g^ + AE (i) x(iT - kT). (1-47) 

The term is found by subtracting the filter response from an ideal 

response d(iT). The factor A is a variable step length which is 
adjusted to improve the filter response in driving toward zero. 

Floating Point Digital Filters 

A floating point digital filter is one which is implemented by 
a computing device which executes floating point arithmetic in calcu- 
lating the filter's difference equations. Both the filter's coefficients 
and the signal variables are represented in the following format 
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x(iT-T) 


x(iT-KT) 





2-47 


F x R E (1-48) 

where F is a fraction expressed in radix R and E is the exponent value. 

R is usually 2 in 16-bit minicomputers but is 16 in IBM 360/370 machines. 
Digital filters are not usually implemented in floating point for several 
reasons. Floating point hardware is slower than fixed point and is 
more costly. Perhaps a more important reason is that floating point 
quantization errors in signal variables can cause system instability 
whereas with fixed point arithemtic is guaranteed to be stable if the 
filter coefficients yield stable poles in the z-plane. 

Optimal Digital Filtering 

Optimal digital filters are filters used to minimize some per- 
formance .evaluation criterion set for the discrete filter. In this 
section, three topics will be presented: 1) the concept of optimi- 

zation, 2) the optimal control law, and 3) state estimation. 

Concept of Optimization [21]. A system may be described by n 
first order linear or non-linear differential equations in the 
independent variables x^, x^, "‘‘x^ Any system can be so described 
by the introduction of the appropriate number of variables, henceforth 
referred to as the state variables. The n differential equations are 

x = f^(x,u,t) (1-49) 


Suppose that a function 
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mT 

V(u) = / L(x,u,t) dt (1-50) 

o 

is to be minimized by choosing the forcing functions u(t) or some 
other system parameters. L represents the performance criterion 
together with any terms which penalize or restrict the use of 
forcing signals. The minimum value of V(u) is termed the cost. 

A linear optimal system has the following characteristics: 

(a) linear differential equations 

(b) the performance criterion has a quadratic form in the 
state variables and forcing functions 

(c) unrestricted forcing functions and state variables. 

Any system which does not possess all three characteristics is non- 
linear. 

Consider the linear system described by the following set of 
first order differential equations: 

& = Fx + Gu 

(1-51) 

Y_ = Hx 

Now it is desired to calculate u(t) (given the initial values 
x(0)) such that the cost function V(u) is minimized. 

It is proposed to approximate the system by a discrete time version. 
The time interval is divided into m equal sub-intervals T and the forcing 
function ia is to be held constant during each subinterval. The system 
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is considered to be described by a sequence of transitions from the 
(k-l)th to the (k)th state. 

Solving the set of first order differential equations, we find 
the following transition equation: 

x(kT + T) = <f>(kT + T, kT)x(kT) + T(kT + T, kT)u(kT) (1-52) 

where 

kx+T 

r(kT + T, kT) - / 4>(kT + T, T)G(T) dT. 

kT 

and $(t,T) is found as follows: 

1. When F is time varying, 3> is computed from 

d [*(t>.T.U = F(t) *(t,T). 
dt 

2. When F is constant $(t,T) = 4>(t-T) is computed by 

t(t-T) = e F(t - T > - l 1F(t ' t ° ) - 1 ! • 

k=o k! 

The relationship between continuous and discrete systems is shown in 
Fig. 16. 

In addition, the integral to be minimized is replaced by the 


summation 




Fig. 16. The Continuous and Discrete Systems 
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m-1 

V(mT) = T [ L[u(kT), x( kT ) > kT l (1-53) 

k=o 

The minimization of V(mT) for discrete systems will be considered 
for two cases: a) the optimal control law, and b) state estimation. 

Optimal Control Law [22], Consider the continuous system equations 
to be of the form, 

x(t) = F(t) x(t) + G(t)u(t) 

X(t) = H(t)x(t) (1-54) 

mT 

V(mT) = x T (mT)Ax(mT) + / x T (t)B x(t) dt 

o 

mT 

+ / u T (t) C(t)u(t) dt 

o 

where A = terminal state weighting matrix 
B(t) = state weighting matrix 
C(t) = control cost matrix. 

The optimal controller is obtained by solving the nonhomogenous matrix 
Riccati equation 

dS Tit 

— = - SF - F S + SGC -1 (0)C i S - B(0) . (1-55) 

dt 

If F and G are constant, 

S (mT) = [<t> 21 (mT) + <J> 22 (mT)A] [<^ u (mT) + ^^(mTjA]" 1 


(1-56) 
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Where 


MmT 
e A 


*11 *12 
*21 *22 


and 


M = 


-F 

B(0) 


-1 T 
GC G x 


Once S(mT) is known, the optimal control vector can be obtained from 


u opt (t) = D(mT - t)x(t) 


(1-57) 


where 


D(mT - t) = C _1 ( t ) G T S (mT - t) . 

In block diagram form, the optimum controller can be depicted as in 
Fig. 17. Notice that to find u^j. the state vector x(t) is necessary 
for calculation of the optimal input. In most systems x.(t) is not 
available; £(t) is available instead. Hence, we "estimate" x(t) using 
£(t) as shown in Fig. 17. 
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State Estimation [23] . It is desired to find an optimal estimate 

A 

x for the state variables x for a system defined in Fig. 17. The 
system output x is measured every T seconds; call the measurement 

z/nT) = x( nT ) + v(nT) 

= H(nT) x(nT) + v(nT) . 

simplifying the notation 


H jc + v 
n— n — n 


where v^ is measurement noise and 

E[v v T ] = R{ 

-n-m t> mn 

E[v ] = 0. 

-n 


(1-58) 


(1-59) 


The estimation scheme is to predict the present value of the state 
vector by using the last predicted value and updating it with the 
present measurement. 

5 n (+) = S„(-> + K n tz n -H n x n (-)] (1-60) 

w h ere x n (+) and x n (-) are estimates of the state vector ^ after 
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and before the measurement s n » at time nT. The is the optimum 
weighting matrix. Let the error in the estimate be 


?n( + > “ - * n 

^ A 

Xn(-) = Xh(-) - x n 

Substituting (1-58) and (1-61) into (1-60) 


(1-61) 


5n(+> " <1 " W5n(“) + Vn * (1 “ 62) 

Define 

p n (+) - E[x n (+) xj(+)] (1-63) 


However, 


E [x n (-)v^] = E[v n x^(-)] = 0 

because of uncorrelated measurement errors. Thus, (1-63) becomes 

V + > * d - - W T + Wn ' (1 - 64 > 

The cost function to be minimized in state estimation is the sum of 


the diagonal elements of the error covariance matrix P (+) : 
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m-1 

V(mT) = l (E[x*(+)x (+)] } 

— n — n 

n=o 

The V(mT) is minimized by K n> The solution is 

K „ ■ p n(-) H il H n p n<->< + ‘‘n]' 1 = V+^X 1 (1 - 65) 

Substituting into P n (+) results in 

P n <+> - P n (-) - P„(-)H P [H n P n (-)Hj + K,,]- 1 V.W, (1-66) 


The equation set for the state transitions of the discrete system 


x 

~n+l 


*n*n + -n' 


(1-67) 


Again, using (1-67) and (1-61) in (1-63) 


w-> 


*n P i +) *2 + Q n* 


(1-68) 


where 


w = r[kT + T, kT]u(kT) 
— n — 

E[ww^] = 06 
-n-m n mn 
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E[w n ] = 0. 

The results given above are now summarized in Fig. 18. 

Nonlinear Filtering [24] 

Reference [24] presents a class of nonlinear systems which obey a 
principle of superposition. In particular, the synthesis of nonlinear 
filters for signals which can be expressed as a product or convolution 
of components is examined. Practical applications in speech and image 
processing are illustrated. 

Range Adaptive Digital Filtering [25] 

In many applications the digital filter's input signal tends to 
dwell near zero with occasional perturbations away from null. Range 
adaptive digital filtering has automatic scaling of its input, internal, 
and output signals to prevent arithmetic overflow. This is a hardware 
concept and will be further examined in PART 3, Mechanization of Digital 
Filters. 

Random Sample Skipping [26] 

In certain time-shared applications of digital filter hardware 
several input/output sequences, say n, of numbers can be handled by 
a single special-purpose computer which looks like n digital filters. 

If the sampling rates for each filter is different, then inevitably 
conflicts for the arithmetic unit will take place and certain input 
samples will essentially be lost. This process can be described as 



Fig. 18. The Continuous System with Optimal State Estimation 
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nonlinear random sample omission. Reference [26] shows that in some 
cases random sample omissions in a closed-loop system with random 
inputs can be beneficial in reducing the mean-square value of a nulling 
error signal. 

Block- Floating-Point Filters f271 

Block-floating-point is a compromise between fixed-point and floating- 
point arithmetic. In fixed point arithmetic no scaling is used for 
addition or multiplication of numbers. In floating-point, automatic 
scaling is performed for each product or sum calculated. In block- 
floating point arithmetic, numbers are expressed as a fraction and 
exponent (as in floating point); however, scaling is performed once 
for an entire expression instead of for each operation. 

For example, 

^ - x n + a^n-i + “ ’ + a N y n -N (1-69) 

would be calculated as 



n 


y n - Vn + a lVln + • • • + a N A n w Nn 

A 

x n = Vl x n 
w in = Vl y n-i 


(1-70) 
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The scaling factors Ajj and A R are powers of 2 and are determined as 
follows 



A" = A nVl 


where Cn is the maximum characteristic of the variables x n , wi^**’, 
w Nn* , the calculations for y n , and w^ n involve only 

A 

scaling (shifting). Once scaling is performed in y n , then all the 
arithmetic calculations are performed in fixed-point. The block-floating- 
point realization is summarized in Fig. 19. 

Sample-Rate Reduction Digital Filters [28] 

There exists a direct relation between input sampling frequency 
and the computational rate of the digital filter hardware implementation. 
In order to prevent input frequency aliasing, two common practices 
are to sample at a high rate or to use an analog low-pass filter before 
the A/D converter. Reference [28] suggests sampling at a high rate 
and using a digital low-pass filter whose output can be sampled at 
a much lower rate to furnish the input signal for some digital signal 
processing system. Advantages include the elimination of phase 
distortions which are inevitable in analog aliasing filters. 








II. TRANSFER FUNCTION SYNTHESIS 


The synthesis of transfer functions for digital filters is 
surveyed in this section. The survey is subdivided into nonrecursive 
filters, recursive filters, and sample designs. 

Nonrecursive Filters 

The synthesis of nonrecursive digital filters consists of determining 
the coefficients hj. of the expression 

M-l _ ± 

H(z) = l h ± z . (2-1) 

i=0 

In the z-plane this amounts to placing zeroes anywhere in the plane 
with all poles falling at the origin. 

Specification of Frequency-Domain Zeroes [29] 

The design of nonrecursive digital filters in the frequency domain 
consists of specifying a finite trigonometric polynomial which satisfies 
some criteria. Here the polynomial is defined by placing its zeroes. 

The frequency charactrristic of (2-1) is defined as 

M-l 

H(f ) = l h e _ J 2lTfn . 
n=0 
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Replacing f by the complex variable <J> = f + ja 

H(4>) = K x n [1 - e ”J 2 ir (W n )] ( 2 - 2 ) 

n=l 

where 

<♦„> - {£„ + j» n ) 

The {<j> n } are the zeroes of H (<(>) in the central period. Hence, the 
scale factor and the M-l central zeroes completely specify H(f). 

The function is factored into stopband and passband zeroes 

H(f ) = Hg (f )H p (f ) (2-3) 


where 


N S 

Hg (f ) = Kg n [1 - e”^ 2Tr(f ”^ n ) ^ 

n=l 

N p 

Hp(f) - Kp n [1 - e~-^ 27 r ^ f ~^n^] 


(2-4) 



2-64 


Once the zeroes have been apportioned between the stopband and passband, 
the passband and stopband zeroes are positioned to give a "good" shape 
for H(f). The procedure is demonstrated in [29], 

Frequency Sampling [30,31,32] 

The technique of frequency sampling may be used to synthesize 
nonrecursive filters as follows: 

1. Choose a set of frequencies at which the sampled frequency 
response is specified. 

2. Obtain the values of the continuous frequency response of the 
resulting filter as a function of the filter parameters 
(defined below) using the sampling theorem. 

3. Compare the interpolated frequency response with the desired 
filter and search for a minimum of some filter characteristic. 

4. When the minimum is found, the parameters are used to realize 
the nonrecursive filter. 

The frequency samples in step 1 are specified in the passband 
and stopband; however, in the transition region several samples are 
left adjustable and these are the parameters used in steps 2 and 3. 

The references [30,31] describe computer aided design programs which 
essentially automate the optimizing process. 

Windowing [19.31.331 

Windowing filters, discussed in Chapter 1, are nonrecursive 
filters whose finite impulse response is found by terminating an 
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= filter coefficients 
P = even integer 

f^ = normalized sampled frequency 

|f k | ^ 1/2 

2W = sampling rate. 

Equation (2-5) may be minimized using the simplex method of linear 
programming. Digital filters designed in this manner are sometimes 
said to have minimax responses. 
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Recursive Filters 

The recursive digital filter has the form 

? _i 

L a i z 

H(z) = i= °— (2-6) 

1 + Zb ,.- 1 

i=l 

The synthesis of recursive filters is the task of choosing the 
coefficients a.^ and b, in order to force the filter to behave in some 
specified manner. 

Direct Synthesis in the Frequency Domain [36] 

The frequency response for (2-6) is found by substituting z = e j2irfT 

l ai e^ 2irfiT 

H(f ) = — (2-7) 

1 + l b ie -J 27rfiT 

i=l 

If N(f) is the numerator of (2-7), then 


|N(f) | 2 - ( l a. 


ie - j2lrfiT )( 


j2frfiT. 
a,e ) 


n 

» H c + 2 l H. cos (2irkTf ) 
k-1 


( 2 - 8 ) 


where 
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n 


H o " l k 


k=0 


\ = I a p a q * 

(p-q)-k P 


Equation (2-8) may be further reduced to 


| N(f ) | 2 = l cos 2k (irTf) , (2-9) 

k=0 

where are constants. Hence equation (2-7) may be written as a 
rational function in cos (irTf) [or sin(irTf)]. Any such rational 
function may be specified by the roots of the two polynomials. 

Consider the Butterworth lowpass filter in the analog domain 


H^f )! 2 = 


1 + (^-) 2P 

fc 


Since the term sin(irfT) corresponds to f in the discrete case 


! h 2 (f ) 1 2 - 


i 


i + 


' sin(TrfT) 1 2p 
sin(7rf c T) 


(2-10) 


represents a lowpass digital filter. To find the filter coefficients 
solve for the roots of the polynomial. An example design in the 
continuous case is presented later in this chapter. 
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Other filter types (bandpass, highpass, band stop, etc) may be 
designed using this technique. 

Sampled Data Transformations [37] 

This section describes a mapping technique for designing recursive 
digital filters. First, a suitable continuous filter G(s) is found, 
and then a mapping function from the s-plane to z-plane is employed to 
find the digital equivalent filter D(z). Hence, first we review 
continuous filter design and then employ the sampled-data transformations. 

Continuous Filter Design . The design of continuous filters can 
be accomplished by first designing several low pass filter transfer 
functions G(s), called prototype or normalized designs; the prototypes 
have a critical or break frequency of one radian/sec. The prototype 
is used to realize a filter for a given specification by using the 
frequency transformations listed below: 


Low Pass: 


Band Pass: 


Band Stop: 


s + s/w u 

s 2 + u u a>£ 
s 

8 <v *0 

s _ 

o 

s* + 0) 0)» 
u t 


High Pass : s -> a> u /s 


( 2 - 11 ) 


where 
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u> = upper cutoff 
u 

i0£ = low CUtoff 

Five prototype filters will be discussed in this section: 
Butterworth, Bessel, Transitional, Chebyshev, and Elliptic designs. 

Butterworth : The Butterworth approximation to the ideal low pass 

filter is defined by the squared frequency magnitude function 

| G(to) | 2 = 1/ [1 + (w 2 ) n ] (2-12) 


where n is the order of the filter. The Laplace transfer function is 
given by 

G(s) G(-s) = 1/[1 + (-l) n s 2n ] , 

or 


G(s) 


n 


n 

j=i 


l 

(s + bj) 


where 

bj = _ e iTr [(l/2) + (2j-l)/2n] i = /T 
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Bessel : The Bessel filter approximation for the linear delay function 

e~ TS may be written 

G(s) - — (2-13) 

V 8 > 

where Kq is a constant term and B n (s) are Bessel polynomials. 

B 0 " 1 
B^ “ s + 1 

B„ - (2n-l) *„.! + s 2 B n _ 2 

l/n 

The roots of B n (s) are normalized using the factor (K Q ) 

Transitional : The transitional filter combines roots of the n th order 

Butterworth and normalized Bessel filters according to a transitional 
factor TF. Let 

Tj •» magnitude of transitional pole 
r^j = magnitude of j** 1 Bessel pole 
= angle of j ^ transitional pole 
0^ = angle of Bessel pole 

0£j = angle of Butterworth pole. 
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the poles of the transitional filter are then described by 

r j ■ r u TF 

6j - + IF(e 1:) - 9 2 j) (2-14) 


Chebyshev : Chebyshev filters exhibit better cutoff characteristics 

for lower order filters than do the above designs. Chebyshev type I 
and type II filters are defined by 



2 = 1 

1 + e 2 X 2 (o>) 
n 


and 


G2(oj) - 


1 + e' 


w 

T n (o3 r /a)) 


1 2 


where 


Tn 


(w) 


cos(n cos"^(u) 
cosh (n cosh~l<i)) 


T 0 " 1 


0 £ u £ 1 
0 ) > 1 


T 1 = 


(2-15) 


(2-16) 


0) 
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T 2 = 2u 2 - 1 
T-j(w) = 4o>^ - 3u> 

The order of the filter n is determined by specifying inband ripple E 
and the lowest frequency at which a loss of a db is achieved. Hence, 

e = ( 10 E / 10 - l ) 1 ' 2 

(2-17) 

A 2 - 10 a/1 ° 

and 


n = COSh ' 1/ A 2 - 1/e > 

cosh”^ ( w ) 

In equation (2—17), the variables E, a, or must be adjusted so the 
n will be an integer. The type I filter differs from the type II in 
that the type I exhibits equiripple in the pass band while type II 
has equiripple in the stop band. 

Elli£tic: The Elliptic filter has equiripple in both the pass and 

stop bands. Hence, this type design usually achieves the desired 
frequency response with a lower order n than any of the above types. 
The elliptic filter is determined by 
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|G(u>) | 2 = 

1 + e 2 i|> 2 (io) 
n 


where 


n odd 




sn[n ^M sn_1 (w ;k );ki] , 

K(k 1 ) . 

sn[K(k 1 )+ N— -±- sn _1 (a);k) ;ki] , 
L K(k) A 


n even 


with 


X 


0) dto 

/ = Elliptic integral of the 

0 [(l-u) 2 )(l-k 2 w 2 )] 1/2 first kind 


sn[x;k] = u) = Jacobian Elliptic function 


K(k) = complete Elliptic integral of the first kind 
tt/2 


o (1 - k 2 sin 2 <j))^^ 2 


k = 1/u) 


k x = e(A 2 - 1) _1/2 


e = (10 E/1 ° - 1) 1/2 


A 2 = 10 a/1 ° 


where e, a, w r were defined for the Chebyshev filter, the order n 


(2-18) 


is found by 
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K(k')K(k) 

KtkpKtk') 


with 


k' - (1 - k 2 ) 1/2 


k' = a - k 2 > 1/2 


The result of any of the five design methods results in a Laplace 
transfer function G(s) for the desired frequency response. 

Sampled- Data Transformations . Once the continuous transfer function 
G(s) has been determined, the transformation to the discrete or z-plane 
is made. Three methods of transformation will be presented: the 

standard z-transform, the bilinear z-transform, and the matched 
z-transform. 

Standard z-transform: The problem of converting a continuous filter 

to a discrete one was presented earlier. It was shown that 


V< s) 

-2 G(z) - Z[G(s) ] 

E ± *(s) 


and that 
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E q (s) 

0,(8) Z[G(s) ] . 

E i *(s) 


But for small T 


G ho (s> 5 1 


Hence , 


E(s) 

-2 = TZ[G(s) ] . 

Ei*(s) 


Define the digital filter D(z) equivalent to G(s) to be 


D(z) = TZ [G(s) ] , 


where Z[G(s)] is the standard z-transform of G(s). Hence, 


D(z) = T l 


\ 


k -! 1 - e -Tb k -1 


(2-19) 


Note that the standard z-transform can be used only on bandlimited 
signals (f < fs/2). 
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Bilinear z-transform : The bilinear z-transform may be used to obtain 

a discrete equivalent of G(s) as follows: 

D(z) = G'(s) 

s = (2/T) (1 - z _1 )(l + z" 1 )" 1 (2-20) 

Where G (s) is a continuous filter whose critical frequencies differ 
from G(s) by 

f' = 1 /ttT tan (irf,,!). (2-21) 

Relation (2-21) is used before the continuous filter G(s) is designed. 

The new filter G (s) is designed instead and then transformed to the 
z-plane by (2-20). The bilinear z-transform is a bandlimiting trans- 
formation with relatively flat magnitude characteristics in the pass 
and stop bands. However, the time response will be considerably 
different. 

Matche d z-transform : The matched z-transform matches the poles and 

zeroes of the discrete function to those of the continuous one. The 
digital equivalent of the G(s) function is calculated as follows: 

D(z) = G(s) 

s + a^ = 1 — z - -*-e _a ^T 


s + bj = 1 - z"^e~kjT 


( 2 - 22 ) 
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If G(s) has no zeroes, it is sometimes necessary to multiply (2-22) by 
(1 + z'V, N is an integer. 

Summary ; The standard z-transform is suitable for only bandlimited 
functions, while the bilinear and matched z-transforms are suitable 
for all filter types. Thie matched z-transform requires G(s) in 
factored form; standard, in partial fraction form; and bilinear, 
in prewarped frequency form. The standard z-transform preserves the 
shape of the impulse-time response; the matched, the shape of the 
frequency response; and bilinear, the flat magnitude gain-frequency 
response characteristics. An example filter is designed and discretized 
in the following example. 

Design Example . In this section a digital filter will be designed 
using the techniques summarized above. 

Suppose it is desired to design a bandstop filter G^(s) with 

w = 200 = 2 tt( 31.831) 
oj£ = 170 = 2ir(27. 056) . 

Multiplied times this filter will be a low pass filter G2(s) with 

= 600 = 2ir(95.493), with a d. c. gain of 1.356. The band stop filter 

will be designed from Butterworth, Bessel, and Chebyshev I prototypes 

with n = 2. The low pass filter will be designed with n = 1. The 

prototype of G 0 (s) = — i — , 

1 s + 1 
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The prototype filters for G^(s) are found below. 


Butterworth : The Butterworth filter is defined by 


Gi (a)) 


1 + of 


G^(s)G^(-s) - — — 


1 + s" 


G, (s) 


1 (. - . .1W4, 


GiU) - L 


+ /2 s + 1 


Bessel ; The Bessel prototype is defined by 


K„ 


G(s) = 


Bo(s) s + 3s + 3 


Chebyshev I : The Chebyshev I filter is defined by 

I G 1 (to) | 2 = 

1 + e 2 T 2 (u)) 

T 2 (co) = 2uj 2 - 1 
e = (10 E/1 ° - 1) 1/2 
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A 2 - 10 a/1 ° 


„ , cosh-V A 2 _ ~ ) 

cosh 1 (oj ) 
r 


Let E = 1.33 db, then 


e = (10* 133 - 1) 1/2 = .5 


Let the filter gain be down 6 db at 


a = 6 


A 2 = 10* 6 = 4. 


cosh ' 1(> 7TT^ > 

cosh 1 (co ) 
r 


n 


= 2 


cosh = -j cosh”^*^!) = .44 


u> r = 1.098 
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Hence, 


|G 1 (oj) 


4 2 

OJ - 10 


+ 1.25 


G (s) = 

(s + 1.057 /31.75 e ) (s + 1.057 /-31.75 0 

G.(s) = -= 

s + 1.308s + 1.118 


The analog filters are designed from the prototypes by setting 


G(s) = G^s) 


X G 2 (s) 


s 


S K ~ a i> 


2 

s 


+ 0) Ojp 

u Z 


S = s/(0 , 
n 


and adjusting the d.c. gain to be 1.356. The resulting filter equations 
are given by 
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Butterworth: 


G(s) 


, s 4 +68000s 2 +1.156X10 9 

813*6 

( s 4 +1272 . 8s 3 +68900s 2 +4 . 32 75X10 7 s+1 . 156X10 9 ) (s+600) 


Bessel : 


G(s) = 2440.8 


s 4 +68000s 2 +1 . 156X10 9 

(3s 4 +27000s 3 +2 . 049X10 5 s 2 +9 . 18X10 7 s+3 . 465X10 9 ) (s+600) 


Chebyshev I : 


G(s) = 909.6 


s 4 +68000s 2 +1 . 56X10 9 

(1. 118s 4 +1177 . 2s 3 +76924s 2 +4 . 0025X10 7 s+1 . 2924X10 9 ) (s+600) 


The filter equations above were plotted for db and phase, <f> , versus 
frequency as shown in graphs 1, 2, and 3. Since the plots are nearly 
identical, the Butterworth G(s) was chosen to be discretized by the 
standard, bilinear, and matched z-transforms, with T = .001. 

The Butterworth design for G(s) may be written in partial fraction 
expantion 

40.567 2.4402X10 -4 + J7.5033X10" 5 

G(s) = — - + 

s + 27.314 x + .35375 + jl85.39 

+ 2.5502X10" 4 - j7.5033X10“ 5 


x + .35375 - jl85. 39 
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1642.1 j -869.03 
+ + 

s + 1244.8 s+600 

The standard z-transform is taken 

a aT 

s + u 1 - e -u ^z - ^ 

and 

a + ib + a - ib T f2a1 + [2e~ uT (bsinvT - acosyT^z " 1 
s+u+iv s+u-iv 1+ [-2e -uT cosvT]z -1 + [e - 2 uT ]z -2 

Hence , 

4. 0567X10 -2 4. 8803X10 -7 - 4.420xX10 _7 z -1 

D(z) « zi + 7 =2 

1 - . 97306z 1 - 1.9654z 1 + .99929z 

1.6421 -.86903 

+ 7 + 7 

1 - . 28800z -1 1 - . 54881z 

The frequency response of this function is found by letting z = e^^. 
The plot is shown in graph 4. Note that this response is entirely 
inadequate. The standard z-transform is accurate only when G(s) 
is limited to frequencies less than 1/2T, or in this case, 500 Hz. 
This condition is violated as is seen in the plot of the continuous 


Butterworth design G(s). 
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The bilinear z-transfora requires a prewarped frequency scale for 
the Butterworth G(s) design, so 

w s = f tan ( ~f T). 



unwarped 

warped 

“u 

200 

200.67 


170 

170.41 

U) 

n 

600 

618.67 


The Butterworth design to be used in this case is 


G(s) = ± 

s 2 + yfl s + 1 


x 1 

s + 1 


s (200.67 - 170.41) 
s 2 + (200.67)(170.41) 


s = s/618.67 


G(s) = 838.92 


s 4 +68392s 2 +l . 1694X10 9 

s 4 +1294.8s 3 +69308s 2 +4.4279X10 7 s+1.1694X10 y )(s+618.67) 


G(s) = 838.92 


(s 2 +3 . 4199X10 4 ) 2 

(s+26. 987) (s 2 +. 70708s+34199) (s+1267 . 1) (s+618. 67) 
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The bilinear z-transform is found by letting 


D(z) = G(s) 



1 - z 
1 + z 


-1 

^1 


d(z) = .19509 (l-i.geeiz^+z'Va+z" 1 ) 

(1-. 97337z 1 ) (l-^9654z _1 +. 99930z -2 ) (1-. 22433z _1 ) (1-. 52749z _1 ) 


The frequency response for this function is found with z = e JwT and is 
plotted in graph 5. Note that this plot closely matches graph 1. 

The Butterworth G(s) may be factored as follows : 


G(s) = 813.6 (s-j!84. 39) (s+1184,39) (s-il84.39) (s+jl84.39) 

(s+27. 314) (s+. 35375+jl84.39) (s+.35375-jl84.39) (s+1244. 8) (s+600) . 


The matched z-transform is given by 


D(z) = G(s) 


s + a = 1 


-aT - 
e 


(s + u + iv) (s + 


u - iv) - 1 - 2e~ u ^cosvTz - ^ + e -2u T z -2 
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and D(l) is set equal to 1.356, the d.c. gain. Hence, 


D(z) 


.34607 


. (l-1.9661z~ 1 +z~ 2 ) 2 

(1-. 97036z _1 ) (1-1 . 9654z -1 +. 99929z“ 2 ) (1-. 28800z -1 ) (1-. 54881z _1 ) 


The frequency response of this function is plotted in graph 6. Note 
that the matched z-transform (like the bilinear) gives a good approxi- 
mation to the response of graph 1. 

Digital Compensators [381 

Digital filters are often employed as compensators for discrete 
control systems. Two common techniques for designing these compensators 
are root locus and Bode plots. 

Root locus . A typical discrete-time closed-loop control system 
is demonstrated in Fig. 20a. Let 
n 

I a i z_i 

D(z) = K izP (2-23) 

I biz' 1 

i=0 


where K is a variable constant and a Q = b Q = 1. The root locus 
technique is outlined below 

1. Find the characteristic equation 


1 + D(z)Z[G ho (s)G(s)H(s)] 



3 CYCLE'S X 1 40 DIVISIONS 


















2-93 


2. Place the poles and zeroes of D(z) inside the unit circle 
in order to make the roots of the characteristic equation 
stable for some range of K. 

3. Vary K from 0 to 00 and solve for the closed-loop roots of 
the characteristic equation. 

4. Choose an appropriate value for K. 

In practice steps 2 and 3 are repeated on a trial and error basis. Once 
the procedure is complete, D(z) in (2-23) is completely specified. 

Bode plots . Bode plots are amplitude and phase plots for a 
transfer function constructed using the asymptotic behavior of simple 
first and second order factors in the numerator and denominator of 
the function D(s). The plots are 

db = 20 log|D(j2irf) | 

♦ = /D(,12irf ) . 

Once the proper frequency response has been found, D(s) may be mapped 
to the z-domain using the bilinear z-transform. 

Frequency Sampling [3] 

Earlier in this report the technique of implementing a finite 
duration impulse response filter in a recursive manner was presented. 

The coefficients must be integer powers of the first one for this 
technique to be applicable. 
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Nonlinear Programming [34] 

Nonlinear programming can be used to design both recursive and 
nonrecursive digital filters. The filter is written as 


H(z) = g 


— 1 

s 1 + a^z 4- b^z 

n lo 

1 + c^z -1 + d^z~^ 


(2-24) 


or 


H(2) - g + l *1 * " . t 2 ' 1 

= ^ 1 M r. ^”1 -i. J 


An error function is formed 


E k = |H(e j27rf k)| 2 - | F(2Wf k ) | 2 k = 1,N 


(2-25) 


(2-26) 


where f^ are the discrete frequencies, 2W is the sampling rate, and F 
is the desired continuous frequency response. Note that E^ is every 
where a differentiable function of a^, b^, c^, d^, and g. The errors 
must satisfy 

- L k± E k± u k k = 1 ,N. (2-27) 


From (2-27) we may define 
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<\ ' » u k - E k 
A 

H k = QL k + E k k = 1,N 


where Q > 0. 

A penalty function such as 


N N 

Q+ l + l ± 

k=l G k k=l G k 


(2-28) 


(2-29) 


is formed. A suitable computer program (such as the Fletcher-Powell 
algorithm [39]) is used to minimize the penalty function with respect 
to Q, g, a^, b^, c^, and d^. Then the factor r is divided by a factor 
and the process is repeated until Q becomes nearly constant. If Q is 
less than unity the procedure stops; otherwise, increase the number 
of stages s of the filter and repeat the above procedure until a Q 
is found less than unity. 

Optimal Digital Equivalent [40] 

In this section the problem of determining an optimal digital 
equivalent D(z) for a continuous filter G(s) is considered (see Fig. 21) . 
The coefficients of D(z) are determined by fitting the input and outputs 
of the two filters. Let 


e od (j) = e o ( j>» 


(2-30) 
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and 


D(z) = 


E ,(z) 
od 

Ej (z) 


a + a, z - ^ +* * *+a ,z' n+1 

o 1 n^l 

1 + b, z -1 +• '*+b z _n 
1 n 


(2-31) 


The problem then becomes one of choosing a £ and b£ in D(z) such that 
(2-30) is satisfied. A difference equation for (2-31) is 


e odt5) l a £ e i(J“^) " l h l e od ( - i ~^‘ (2-32) 

l=o 1=1 

Substituting (2-30) into (2-32) 

n-1 n 

e(j) = l a/eiCj -£) - [e 0 (j) + l b,e 0 (j-£)] (2-33) 

£=o £=1 

where e(j) is driven to zero by minimizing e2(j), the mean squared 
error. 

In vector form (2-33) becomes e(j) * £ T (j)c - e 0 (j) where, 


c T = [ao»-**»a n _ 1 ,-b 1 ,***,-b n ] 

£ T (j) = I e i (j) ’ * ’ »e^(j-n+l) ,e 0 (j-l) , • • • ,e 0 (j-n) ] . 


(2-34) 


The mean square error is 


— N 

e 2 = lim 1/ (2N+1) £ e 2 (j) 

N-*» j=~N 


(2-35) 
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Equation (2-35) is minimized by 



= 2 lim 
N-*x> 


N 

1/ (2N+1) l a(j)e(j) = 0 
j=-N 


or 



N 

N 

(lim 1/2N+1 

l a.(j)q'(j))c = lim 1/2N+1 

l i(j)e 0 (j) 

N-KO 

j=-N N-x» 

-1 1 

j=-N 

1 


R 

r 


and 


c = R -1 r 


(2-36) 
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Equation (2-37) and (2-38) may be written 


A B 

R = [ ] r T = [E F]. (2-39) 

C D 


The elements of (2-39) are of the form 


A: $e^ei(kT) 

B: $ ei e 0 (kT) 

C: $e oei (kT) 

D: 4>e 0 e 0 (kT) 

E: $eie 0 (kT) 

F: $e 0 e 0 (kT) (2-40) 


where 


N 

*xy(kT) - lim 1/(2N+1) £ x(j)y(j-k); 

N-*» j=-N 

the $ X y is the correlation function for discrete sequences. Since the 
input signal power spectrum <J»e 1 e i (s)and the analog filter G(s) are 
known, (2-40) is determined by 


4 > ei e i(kT) * $ e l e l(T) 


-1 


x = kT 


[$e i e i (s)] 


t = kT 


-1 


4> e i e 0 (kT) = [3>e i e i (s)G(s) ] 


t = kT 


$e 0 e i (kT) , 


(2-41) 


and 
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* e o e o (kT) 


1 [<J>e i e i (s)G(s)G(-8) ] 


t = kT. 


The digital filter determined above should have higher order than 
its analog counterpart so that the mean squared error will be small. 

Sample Designs 

In this section some example digital filters are listed. 

Bandstop Filter 

A digital bandstop filter was designed earlier in this chapter: 


D(z) = 


. 34607 (1-1 . 9661z -1 +z -2 ) 2 


— T" 

(1-. 97036z _1 ) (1-1 . 9654z _1 +. 99929z" 2 ) (1-. 2880z _1 ) (1-. 54881z~ ) 


(2-42) 


The frequency response for T = 0.001 is shown in graph 6. 

Digital Resonators 

A digital oscillator is formed by placing complex poles on the 
unit circle: 

D(z) = (2-43) 

1-2 cos(2 fT) z -1 + z -2 


where T is the sampling period and f is the frequency of oscillation. 
Experimental results are available in [41]. 
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Digital Differentiators [42] 

The differentiator is a necessary part of many practical systems. 
The digital differentiator may take many forms; perhaps the best is 
a forth order recursive design shown below: 


_ A (1 - az 1 )(1 - bz ^(l - cz ~ dz 1 ) 

D(z) - A — — 

(1 - ez”^) (1 - fz _1 )(l - gz _1 )(l - hz’ 1 ) 


where 


A = 0.36804011 
a = 0.99999949 
b = -0.86810806 
c = 0.32672838 
d = -.44183252 


e - -.10779165 
f = -.87602073 
g = 0.33494085 
h « 0.51312758 . 


This differentiator was designed using nonlinear programming. 

A nonrecursive wideband differentiator can be constructed for N 
samples by the relation 


^ = k/(N/2) 


k = 0, N/2 


= (N-k) / (N/2) 


k = N/2 + 1, N - 1. 


(2-45) 


If the center samples are adjusted to optimally minimize the magnitude 
error for N = 16, then 



Gjj/2 = 0.92890015 

G (N/2)-l = 0.86994255 (2-46) 

G (N/2)-2 " °* 75000000 

yields a peak error magnitude of 7 x 10”^ for an 80% bandwidth. 

Low-Pass Filters [43-45] 

Reference [43] presents some 9 example nonre cursive low-pass filters 
of order 11. The designs are found using prolate spheroidal functions, 
least mean-square error, Fourier coefficients, windowing, binary 
weighting, and minimax techniques. The reader is referred to Table 1 
of [43] for the appropriate coefficients. 



III. COEFFICIENT QUANTIZATION 


General 146] 

One effect of finite wordlengths in digital computers is that the 
filter's parameters, or coefficients, must be chosen from a finite set 
of allowable values. Classical design procedures yield filter transfer 
functions with coefficients of arbitrary precision which must be altered 
for implementation using digital computing devices. One approach to 
this problem is to select a filter structure (programming form for 
the difference, equations) which is not sensitive to coefficient 
inaccuracies. For example, realizing a filter directly allows a 
greater chance for instability than cascading or paralleling second 
order modules because it is well known that the roots of polynomials 
become more sensitive to parameter changes as the order of the polynomial 
increases. 

Any programming form, or structure, produces a grid of allowable 
pole/zero locations in the z-plane. The proper structure to choose is 
one for which the grid is most dense in the areas at which the poles/ 
zeroes must be placed for a particular design. It is obvious that 
arbitarily rounding or truncating denominator coefficients could cause 
poles to migrate outside the unit circle causing filter instability. 

Instability Thresholds 147] 

For low-pass filters , a measure of the number m^ of bits required 

to represent the coefficients of a stable filter may be expressed as 

2-104 
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m b = i - N log 2 (2TT BT) 


(3-1) 


where 

B = minimum attainable bandwidth 



for the direct programming form. For the cascade form 

2 - ( m ,, - 2)/2 


(3-2) 


These stability thresholds are valid for filters designed using direct 
synthesis in the frequency domain for sine and tangent Butterworth low 
pass filters. The results may be extended to other filter types. 


Reduced Coefficient Wordlengths [48] 

The cost of implementing a digital filter via a special-purpose 
computer is directly related to the wordlength of its coefficients. 
However, a short wordlength can cause large deviations in pole/zero 
placement. Hence a compromise must be found. The following procedure 
represents one solution to the problem. 

Let the transfer function of the digital filter be 


H(z) 


m 


l v 

i=0 1 


-i 


m 


l c i z 

1=0 1 


-i 


(3-3) 
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where c =1. If we examine the desired frequency response H around 
o w 

the unit circle 



jwT 


)| 


1 in passband 
0 in stopband 


(3-4) 


and is unspecified in the transition regions. If the maximum passband 

and stopband deviations are defined as 6 and 6 

P s 

|H w (e Ja,T )| - |H(e ja,T )| < 6(e ja>T ) 


or 


e(e ja)T ) = 6 1 1 - |H (e ja,T )| | in PB 

P n 

I , ja»T . (3-5) 

= 6 g |H n (e J ) | in SB 

where t is the normalized error function and H (z) = K H(z), a normalized 

n ’ 

transfer function. 

The design of H(z) minimizes max e such that 

max e < 1 (3-6) 

using standard minimax procedures. If (3-6) holds for a set of parameters 
a, then justification for searching for a second set a/ of reduced word- 
length which also satisfies (3-6), where 
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The coefficients are usually found for the cascade or parallel form. 

The search for a new set a' follows a modified univariate procedure 
which is described below: 

1. Several sets of parameters, say 10, are stored in order of 
minimum max c. 

2. Perform a univariate search on the best set a^. If no improvement 
can be found, try ay. 

3. Stop the procedure when no better improvement is found for 
any stored coefficient a^. 

Generally, rounding of the coefficients is first performed. A 
univariate search reduces max e by 25 to 50% over rounding, while a 
modified univariate search produces the best results reducing max e 
by 25 to 50% over the univariate search. 

In general, the development of synthesis procedures for quantized 
digital filter coefficients remains an active area of research. 



IV. NONLINEARITIES IN FIXED POINT ARITHMETIC 


In digital computer implementations for digital filters, the 
restriction of finite wordlength produces several nonlinear phenomenon. 
Quantization occurs at the input sampler and in the internal arithmetic. 
Saturation and overflow also manifest themselves. Inaccuracies in 
coefficient representation has been discussed previously. Other noteable 
effects which must be examined are limit cycles and deadbands. 

Quantization Errors 

A digital filter specified by equation (3-3) is implemented by pro- 
gramming constant coefficient linear difference equations. The program 
for the difference equations will consist of the arithmetic operations, 
multiplication and addition (subtraction), and data transfer operations. 

The arithmetic unit of the computing device must be furnished binary 
numbers for the coefficients and variables of the difference equations. 
Since each coefficient and variable is represented by a finite number of 
binary digits, the binary numbers supplied to the arithmetic unit are 
quantized versions of the real numbers expected in the difference equation. 
Hence the digital filter introduces quantization errors into the system 
of which it is a part. 

Quantizer Types [49] 

Signal amplitude quantization results from A/D conversion of the 
digital filter input signal, and from arithmetic operations with in the 
computing device itself. Three common types of arithmetic quantizers 
are shown in Fig. 22; the step-length of each quantizer is h. Fig 22a 
illustrates the quantizing characteristic for a roundoff quantizer. The 
roundoff quantizer approximates the input signal e^ by the closest quan- 


tized value e.j^ as follows: 
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for e^ > 0 
for e^ < 0. 


(4-1) 


Therefore, the maximum error magnitude is ^ . The properties of the 
truncation quantizer is shown in Fig. 22b. This quantizer is less diffi- 
cult to implement than the roundoff type; however, the approximation e iq 
is less accurate: 


0 <_ e^ - e^ q < h for e^ > 0 
-h < e^ - e^q <_ 0 for e^ < 0. 


(4-2) 


Here the maximum error magnitude is h. 

The third arithmetic quantizer presented in Fig. 22c is labeled LSB-1. 
In LSB-1 the least significant bit of quantized binary words is always 
set to "one." For this quantizer e^ q is never equal zero. 

-h < e 4 - e in < h for e., > 0 

“ 1 q 1 (4-3) 

-h < e^ - e£ q <_ h for e^ < 0. 


Again the maximum error magnitude is h. 

Signal amplitude quantization at the A/D converter usually takes two 
forms. If the A/D converts the input signal magnitude to binary form, 
then the quantizer characteristic of Fig. 22b for truncation adquately 
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describes the effect of the A/D. However, if a bipolar A/D is used, the 
bipolar property is usually obtained by an offset bias voltage which 
causes the bipolar A/D quantization characteristic shown in Fig. 23. For 
this quantizer 

0 1 e i " e iq < h (4-4) 

and the maximum error is h. 

In summary, the maximum error magnitude introduced by a quantizer 
at a sampling instant is 

Roundoff: h^ = h/2 

(4-5) 

Others : h^ » h 


The quantizers of Figs. 22 and 23 may be represented in a system as 
an additional input error signal; this process is shown in Fig. 24. Using 
this model for the quantizers, their effect on system response will now 
be considered. Mathematical analysis of quantizing errors may generally 
be described as steady-state analysis, statistical analysis, and error 
bound analysis. Each of these analysis techniques will now be presented. 

Steady-State Analysis [50] 

The steady-state analysis may be divided into three steps. First, 
find the z-plane transfer functions T^ (z) from the jth quantizing error 
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Fig. 23. Bipolar A/D. 


(A) THE QUANTIZER 


e ± (kT) 


Fig. 24. 



(B) MATHEMATICAL REPRESENTATION 


Mathematical Model for a Quantizer. 
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source Nj to the system output, e Q . The total number of quantizers in 
the system is s. Hence, 


E on (z) = l T j< z)N 1 (z) ’ 
j=l J 


(4-6) 


where E on (z) represents the output due to quantization errors. 

Second, assume each error source is a step input of the maximum 
error amplitude h^ for the type quantizer being analyzed. Therefore, 


Nj(z> - 


5 


1 - z 


-1 


(4-7) 


Substituting (4-7) into (4-6) 


s T (*) h! 

w*> - I . - — 

— 7 ! 


(4-8) 


j=l 1 - z 


Lastly, apply the z-transform final value theorem [y(«°) = lim(l - z - *) 

z-*l 

Y(z) ] to (4-8); thus, 


e on ( eo ) = 11m l T.(z)h' 

*+l j-1 J 3 

= l J lim T^ (z)1 h' 

j-1 Lz^i J 3 


If one defines 


K 


ssj 


lim T. (z) 
z-KL 3 


(4-9) 
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then 


e («>) 
on v 7 


J-l 


K ssj h j 


(4-10) 


Equation (4-10) may be used to evaluate the effect of each quantizer on the 
system output under steady-state conditions. 

Another technique for finding the Kg S j weighting constants for (4-10) 
is derived as follows. The standard z-transform for tj (t) is 


T|(z) = l t,(kT)z" k 
J k=0 J 


where kT represents a sample instant. Hence (4-9) becomes 
00 

K . = l t . (kT) 

SSj k=0 j 


(4-11) 


If tj(kT) tends to zero as kT gets large, say NT, 

N 

Kssj “ I fc j< kT > (4-12) 

k=0 


may be used in (4-10) to calculate the steady-state error. The terms 
tj(kT) in (4-12) may be obtained from a simulation of the system by apply- 
ing Nj $z) * 1, a discrete impulse function. Note that the weighting con- 
stants are functions of the system characteristics and not of the quan- 


tizers. 
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Statistical Analysis [51] 

If the input signal to a roundoff quantizer Qj has a dynamic range 
of more than three step intervals hj , the effect of the quantizer may be 
determined by replacing it with a unity gain and an additive white noise 
nj (kT) (see Fig. 24) with a rectangular amplitude distribution density 
function p(nj) of bounds +hj/2 and height 1/hj. The LSB-1 quantizer can 
also be replaced in this manner with p(nj) bounded by +hj and l/2hj . The 
truncation quantizer cannot be represented exactly in this manner, but 
this technique does give a good approximation with p(nj) bounded by 
+hj and l/2hj. Let us continue by analyzing the roundoff case which 

can be easily extended to the others. 

2 

The variance of this rectangular distribution is 


°nj * [_ n jP (n j> dn j “ 


(4-13) 


When the dynamic range of the input signal is greater than three 
quantization levels, the noise input of the quantizer is essentially un- 
correlated between successive sampling intervals, and the autocorrelation 
of the quantization noise becomes 


CO 

Vn>> " 1 ( T - M)/T M < T 

J J 00 J 


(4-14) 


= 0 


t| > T 
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The sampled power density spectrum is defined by 


! *„ „ < nT >* ■ c,n i 

j j n= “°° j j 


12 


(4-15) 


The mean-squared error output due to one quantizer error is 


(kT) = l/2iri J $ n ^(z) Tj (z)Tj (z _1 )dz/z 


"on 


(4-16) 


where T^ (z) « E on ^ z -'/ l ' , j 


(z)/Nj (z) , T is the unit circle, and i = Substi- 

tuting (4-15) into (4-16) and assuming that the total rms output error is 
bounded by the sum of the s rms errors due to the quantizer inputs yields 


t e on^rms — ^ K stj 


(4-17) 


where 


Kstj 


■ f— f 

|^24iTi J 


Tj (z)Tj (l/z)dz/z 


- 1/2 


(4-18) 


The integral in (4-18) may be evaluated by calculating the residues of the 
integrand. 

Another technique for calculating the mean-square output error is by 
using the following identity: 


oo 

-L- f F(z)F(l/z) dz/z = \ f 

2iri J k=0 


(kT)' 


Hence (4-18) becomes 


K 


stj 


l ^ (kT)2 ] 


1/2 
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where t j (kT) Is found from the Impulse response In a simulation of the 
system. If tj(kT) converges to zero for k large, say N, then 


K 


stj 


[t-2 X V-> 2 ] 


1/2 


(4-19) 


This relation may be used instead of (4-18) for many applications. 

Equations (4-18) and (4-19) are for roundoff only. They should be 
altered by substituting hj = 2hj into (4-17) for the general case. 

Quantization Error Bounds [52 1 

Consider an n^ order system described by 


x(k + 1) - Ax(k) + Dr(k) 
e Q (k) = £ T x(k) + d^r(k) 


(4-20) 


where j:(k) is a vector of the system inputs and e Q (k) is the output. 
The sampling interval T has been eliminated for convenience. The 
introduction of quantizers into the system results in 


x^(k + 1) = Ax^(k) + Dr(k) - Bci/k) 

e **(k) ■ c^x^(k) + d^r(k) - f^q.(k) 
o — — 


(4-21) 


where £_(k) represents a vector of the s quantizer error inputs n^ (k) . A 
state variable representation for the quantization error 
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v(k) = x(k) - x q (k) 


e on< k ) = e o (k) ‘ e o q < k > 


(4-22) 


results from subtracting (4-21) from (4-20) 


v(k + 1) = Av(k) + BcL(k) 


e on ( k ) = £ T v(k) + f^g/k). 


(4-23) 


The general solution for (4-23) is 


N ~1 

v(N) = A w v(0) + £ A £ Bcl(N - 1 - £) 

£=o 

N-l 

e on (N) = £^A N v(0) + l £ T A A B£(N - 1 - Z) + f\(N) . 

l=o 


(4-24) 


For N large, 


s / N-l n \ 

“ (N) = j l\ ^ A^jq^N “ 1 " 

e on (N) = l ( I sVkjW 11 - 1 - £) + f\(N). 
j=l\£=o / 


where Qj (N) is the j Quantizer and bj is the j ^ column of the nxs 
matrix B. Since, if a = be + de, la I < lb lx I c 1+1 d lx |e I 


(4-25) 
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vooi i l (Y l A ^jl)lqj (N - 1 - 1) \ 

j-lU«o J / J 


and 




where h! - I q . (n - 1 - l ) I is given in (4-5). Similarly, 
j 1 j max 




'on 


In another form, 


|v(N) < £ mjhj 


l e o„<»>lma*i IW-j 


(4-26) 


(4-27) 


(4-28) 


where 


£j = l |aV|, j - 1, s 

J l=o ~ 



(4-29) 


Note that (4-29) gives weighting vectors mj and weighting constants 
which are functions of the system and not of the quantizers. Hence, (4-28) 
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and (4-29) are useful In helping to choose quantization error schemes for 
systems with digital filters. 

A second method for bounding the output error due to quantizer Qj 
is from the transfer function 


E (z) 

T j (z) = , V 

J Nj (z) 


The impulse response is found with (z) = hj. Therefore, 


E on (z) - Tj (z)hj - [j 3 


(4-30) 


To calculate the worst case output error e Qn due to quantizer Qj 


on 


»Ui[J 0 -|*a« | ] h 5' 


(4-31) 


Similar to the argument employed for equations (4-12) and (4-19) 


l e on (N > L« 1 ^ 


where 

Ku bJ = Jo I c j< k > (• 3 " s - ( 4 - 32 > 

Equation (4-32) may be used to calculate the weighting constants K 
instead of (4-29) . 



A s umma ry of the results of the quantization analysis presented in 


this section is displayed in Table 1. 


TABLE 1: Quantization Analysis 


s 

Figure of Merit » Z Constant., h.j 

j-1 

Analysis 

Method 

Figure of 
Merit 

Constantj 

Steady 

State 

Steady 
State Error 

^ssj = liin T j (z) 
z + 1 

vi : ! *j <ki) 

k=0 

Statistical 

Root mean 
square error 

1 

im 

Error 

Bound 

Maximum 

error 

*ubj - to IcVbjI+lfjl 
•Uj 5 Lo 1 CkT> 1 


Open-Loop vs. Closed-Loop [53] 

The quantization analysis procedures above are equally applicable 
to open-loop or closed-loop systems. However, open-loop analysis of the 
digital filter itself is perhaps the easier approach. It has been 
shown in [49,53] that open-loop analysis can give satisfactory results 
even if the filter is to reside in a closed-loop system. 
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Limit Cycles and Deadband Effects [46,54,55] 

Consider the digital filter 

y n = x n + B y n-l * B " '0.5 (4-33) 

implemented in fixed-point arithmetic with roundoff quantization. If 
the input x n is a impulse function of value 7/8 

y 0 - 7/8 
yi = 1/2 

y2 = 1/4 (4-34) 

y3 = 1/8 

y n = 1/8 n ^ 4 

is the resulting output sequence. Ideally the output should go to 
zero. This type error is called a limit cycle, and the amplitude 
intervals within limit cycles are called deadbands. The deadband for 
(4-33) is 

IjviI - ebn-il i < i >2‘ b 

where b is the number of magnitude bits. Hence 



1 - | 3 


(4-35) 
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For the second-order filter 

y n = x n - Vn-l " Vn-2 (4 ' 36) 


the deadband is 



(4-37) 


The deadband for higher order filters is directly dependent upon 
the programming form. In general, the parallel form yields better 
results because one need not be concerned with the ordering of cascaded 
sections [54], 


Saturation and Overflow [56] 

When a filter is implemented in one’s or two’s complement arithmetic 
and signal values exceed the finite register length upper limit, a 
overflow condition occurs and the results usually changes sign. This 
condition can cause large limit cycles, called overflow oscillations, 
to be excited. These oscillations may be avoided by using saturation 
arithmetic as designed in [57,79]. One must be wary of this solution 
for in many closed-loop control systems, saturating the signals causes 
system instability. Saturation changes the filter output which 
effectively alters, temporarily, the transfer function. 
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Dynamic Range [46 ] 

The dynamic range of a binary signal x n of b + 1 bits is 

0 1 |xj < 2 b - 1. (4-38) 

Increasing the number of bits by one doubles the dynamic range. As 
seen in the last section, it is important that the dynamic range of 
a digital filter in many applications never be exceeded. Hence, several 
techniques may to employed to find b. 

One technique finds the least upper bound on the signal x n and 
uses (4-38) to specify b and hence this limit can never be exceeded. 

More practical solutions use simulation of the filter with typical 
inputs to define the dynamic range of the internal variables. Some- 
times statistical methods are used for non-deterministic input signals. 



V. NONLINEARITIES IN FLOATING POINT ARITHMETIC 

In the past there has been little emphasis placed on research and 
analysis of quantization errors at the output of a floating point filter, 
the reason for this being that most filter implementations use fixed 
point arithmetic. Sandberg [58] was the first to study quantization 
error analysis for floating point filters with [59-62] being more recent. 

As in the case of fixed point filters, quantization error for floating 
point filters has three sources due to finite word length. They are 

1) the quantization of the input signal Xjj into a set of discrete 
levels; 

2) the representation of the coefficients of the filter, a k and 
b^, by a finite number of bits; 

3) the accumulation of roundoff errors caused by arithmetic 
operations. 

Notation 

If we assume the ideal output of the filter is w n and the actual 
output y n , the error at the output of the n^ sample e n may be defined 
as 


e 


n 


y n “' w r 


(5-la) 


where 
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M N 

w n = l a k x n-k - l b k w n-k (5_lb) 

k=0 k=l 

Before the effects of the above error sources are discussed, the repre- 
sentation of floating point numbers with a fixed number of bits should 
be considered. 

A floating point number is written in the form (sgn)2^*a, where b 
is a binary integer called the exponent and a is a fraction between 
1/2 and 1 called the mantissa. As expected, the range of numbers that 
can be represented is determined by the number of bits of the exponent. 

In order to represent a number v in floating point form with a t-bit 
mantissa, the smallest integer exceeding log 2 v is first determined. 

This number is denoted by riog 2 vl . The binary expansion of the fraction 
v/flog 2 vl is then rounded to t bits. If (v) t denotes the t-bit mantissa 
floating point approximation, it is seen that 

(v) t = v (1 + e) (5-2) 

where the error is bounded by -2 -t _< e < 2 t , or [-2,2). 

Error Sources 

Both addition and multiplication in floating point arithmetic 
introduce roundoff error. Let (v^ # V 2 ) t and (v^ + V 2 ) denote, 
respectively, the actual computed product and sum of two numbers V]_ 
and V 2 ; then 


\ 
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( v l* v 2>t = + ( 5_3 ) 

(v-^ + V2> t = (v^ + V2> (1 + e) (5-4) 

where the errors 6 and e are bounded by t-2 -t , 2 -t ). 

The above errors will be regarded as random quantities and they 
will be uniformly distributed in their range (-2 -t , 2 _t ). Making these 
and the above assumptions, a statistical approach will be discussed 
which predicts floating point quantization errors. 

First, consider the effect of input quantization. Supposing 
the quantizer has equal step size h, the input to the filter is 
^ + e^ where each e^ is bounded by -(h/2) <_ e§ < (h/2). Since the 
filter is linear, the output is the sum of the two components, x n 
and e^. In determining the effect of input quantization, e§ is 
considered as white noise with a zero mean and variance h 2 /12. The 
steady-state output component due to e^ is a zero-mean wide-sense- 
stationary (w.s.s.) sequence with power spectral density 

H(z)H(l/z)(h 2 /12) (5-5) 

where H(z) is the transfer function of the filter as repeated below 
M N 

H(z) = ( l a.z~ k ) / (1 + l b,z -k ) 
k=0 k=l k 


(5-6) 
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The effect of coefficient inaccuracy on roundoff accumulation has been 
ignored. 

An expression for the mean-squared value of the error at the 
filter's output due to input quantization is obtained by integrating 
the power spectral density (Equation (5-5)). It is equal to 

l/2irj | > lH.(z)H(l/z)(h 2 /12)]/z dz (5-7) 


Coefficient Quantization 

Considering the effect of coefficient quantization, it is seen 
that each coefficient is replaced by its t-bit representation 
according to (5-2). This means the coefficient a^ is replaced by 
( a k) t » which equals a k (l + ot k ) , with a k bounded in absolute value 
by 2 -t . Likewise, each b k is replaced by (b k ) t which is b k ( 1 + g k ) . 
Because of this, it is abvious that the filter characteristics will 
change. The problem can be approached in several ways. The first, 
and the simpliest, is to compute the frequency response of the actual 
filter with t-bit rounded coefficients by using the actual transfer 
function 

M N 

[H(z)J t = ( l (a k ) t z" k )/(l + l (b k ) t z- k (5-8) 

k=0 k=l 

and then comparing the result with the ideal response of the original 
design. 
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Coefficient rounding can cause movements of the poles and zeroes 
of the transfer function. When this happens, network sensitivity 
theory can be applied to study the changes of the filter response. 

If the poles of H(z) are at z^, i * 1, N, and the poles of [H(z)] t 
are at z^ + Az^, it can be shown that 

N N 

Az i = l I(z£ + L )/ n (1 - Cz i /z n ))]/Aa k (5-9) 

k=l m**l 

m?*i 

where Aa k is the change in the coefficient a k . Likewise, results can 
be obtained for the movement of the zeroes. The change in the filter 
response can be studied from these movements. 

Instability of a filter may occur, due to coefficient error, when 
a filter has poles that are close to the unit circle in the z-plane. 

The problem can be serious when the sampling rate of the filter is 
relatively high, even for low order filters in the direct form. 

Kaiser [62] has demonstrated that for an N^-order low-pass filter 
operating at a sampling rate of 1/T with distinct poles at e - Pk^, 
stability is guaranteed if the number of bits used m^ satisfies the 
inequality 

N 

% > [-log 2 [5yij /2N+ 2 ) ( n p k T)]j (5-10) 

where the bracket denotes the samllest integer exceeding the quantity 
inside. It is also possible to extend the result to include multiple 
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poles and to derive similar results for filters of other than low-pass 
type. 

The effect of coefficient inaccuracy is more pronounced for a 
high-order filter when it is realized in the direct form than when 
it is realized in the parallel or cascade form, which suggests the 
parallel or cascade form should be used for high-order filters when 
possible. Further details on coefficient quantization are given in 
Chapter 3. 


Output Error 

Roundoff accumulation error for floating point filters [59-61] 
is quite different from that of fixed point filters and consequently 
will be treated with more depth than that of fixed point. The errors 
introduced are relative to Equations (5-2) , (5-3) and (5-4) . The 
calculation of the statistical mean-squared error at the output will 
be discussed for the direct programming form with the understanding 
that extension to other forms is easily accomplished [61] . 

It has been shown that for floating point arithmetic the actual 
filter coefficients are a^l + ct^) and bk(l + 6^) where and 6^ are 
bounded in absolute value by 2 _t . The actual computed output y n is 
given by 

M N 

y n = l a kd + a k> x n-k l V 1 + 6 k>y n -k ] 

k=0 k=l 


(5-11) 
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where fZ[ ] denotes the actual computed result by floating point 
arithmetic of the quantity inside the brackets. It is assumed that 
the computation of (5-11) is carried out in the following order: the 

products ajt(l + 80(1 b k(l + &k^ y n-k are * irst formed; the two 

sums are then calculated; and finally the difference is taken to give 
y n . Each of these arithmetic operations introduces a round-off error 
which is characterized by (5-3) and (5-4). A flowgraph of this operation 
may be drawn, as is shown in Fig. 25, which includes all the roundoff 
error introduced in the calculation of y n . From Fig. 25, it is seen 
that <S n> k is introduced when the product of a^l + a k^ x n-k * s f orme d> 
and 6 n ^ is introduced when the computed products of a Q (l + Oq)^ and 
a^(l + otpXjj-.i are added. The actual output y n is then 


M 

■ I “kd + “k> 6 n,k*n-k 


I b k (l + 6k)*„,kyn-k 

k=l 


(5-12) 


where 


M 


®n,o 

(1 + 5„)(1 + S n>0 ) 

n 

(1 



l-i 




M 


®n,k 

a + e n )d + 5 n>k ) 

n 

(1 



i=k 




N 


♦n,l " 

(1 + «„> a + V 1 

n 

(1 



i=2 




N 


^n,k = 

(1 + 5 n>d ♦ * n>k > 

n 

i-k 

(1 


+ ? n,l> 

+ C n>i ), k = 1, 2, * • * ,M 
+ ^n.i) 

+n n ,i), k = 2 , 3, • • * ,L (5-13) 
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The quantities 6 n>k ; 4 n)k ; n n>k ; e n>k ; and are the errors introduced 
at each arithmetic step and they are independent random variables uni- 
formly distributed in 2 _t ) . 

From Equations (5-1) and (5-12) it can be shown that the error 
e n satisfies the following equation: 


X 


= u’ + u' 
n r 


(5-14) 


where b Q = 1 and 

M N 

< ' J 0 Wn-k ' VkVk 

“n * I Wk - l) Vk - ! b k ( *n,k - nVk <5 - 15 > 

k=0 k=l 

In the above equations is due to coefficient rounding; u|J is due 
to roundoff accumulations and the input Xjj is zero mean and w.s.s. 

Both components u^ and uJJ have zero mean and are w.s.s., and they are 
uncorrelated, this being because 0 n ^ k and 4> n#k have a mean equal to 1 
and are independent of ^ and w n . 

From Equation (5-14), the error sequence e n is zero mean and w.s.s. 
with a power spectral density related to those of u^ and u^ by 


$ee( z ) = Il/(D(z)D(l/z))] [4> u * u » (z) +$ u » u "(z)]. 


(5-16) 
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Vu'(z) is calculated from Equation (5-15) and is given by 


* » «(z) = lB(z) - H(z)A(z)]-lB(l/z) - H(l/z)A(l/z)]$ xx (z) (5 


where H(z) (Equation (5-6)) is as previously defined and 


A(z) = l bj k f k 
k=l k 

M 

B(z) = l a k a k z _k (5 

k=0 


Concluding from Equations (5-13) j and (5-15), uJJ is white noise with 
power spectral density as follows; 


$ „ (z) = q 2 /2Trj l(F(z) + G(z)H(z)H(l/z) 

u u J 

-N(l/z)[D(z) - l]H(z) 

-N(z)lD(l/z) - l]H(l/z))$ xx (z)/z dz (5 


M 


n 

where N(z) = l a k z is the numerator of the transfer function in 
k=0 

Equation (5-6) and 


M M 


Hz) - l l Vi F k i 

k=0 i=0 k ,:L 


k-i 


- 17 ) 


-18) 


-19) 
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N N 

G(z) = l l 

k=l i=l 


b k a i G i,i Z 


k-i 


M + 2 - max(k,i) , 
k>1 M + 3 - k. 


k^iork=i=0 
k = i t 0 



N + 2 - max(k,i) , 
N + 3 - k. 


k^iork=i=l 
k = i 4 1 


(5-20) 


The mean squared value of the error e n can now be calculated from 
<J> ee (z) by using 

E{e^} = 1/2tt j £ $ ee (z)/z dz. (5-21) 



VI . PROGRAMMING FORMS 


The structure of a digital filter is described by a unique set 
of constant coefficient linear difference equations. These difference 
equations constitute the digital filter’s programming form. As a general 
rule, for any programming form the lower the order n of the filter 
transfer function the smaller the error introduced into the system by 
coefficient and signal amplitude quantization. Consequently, a n t ^ 1 
order filter is usually factored into second-order modules which are 
paralleled or cascaded to realize the higher orders. The second-order 
is chosen so that complex poles and zeroes are realizable. 

The z-transfer function for any second-order module may be expressed 

a Q + a^z~^ + 

D(z) = (6-1) 

1 + b^z + b 2 

The eleven programming forms presented here will be for the second-order 
module of equation (6-1). For a higher-order digital filter, the following 
procedure applies: 1) Section D(z) into second-order modules, 2) analyze 

each module using the computer-aided design procedure to be developed 
later, and 3) cascade (or parallel) the resulting designs to realize the 
original D(z). 

This section will summarize, for equation (6-1), eleven different pro- 
gramming forms and the attributes of each needed for quantization analysis 


2-137 



2-138 


by steady-state, statistical, and upper-bound techniques. In particular, 

the transfer function Tj (z) from the j th quantizer to the filter output 

for equation (4-9) and (4-18), and the discrete-time difference equations 

til 

necessary for the impulse response from the j quantizer to the filter 
output for equations (4-11) , (4-19) , and (4-32) , will be listed for each 
programming form. The tabulation of the eleven programming forms is 
a result of [38, 63-66]. Many others are possible as seen in [67-70,76]. 

Direct Form 

The direct programming form for equation (6-1) is shown in Fig. 26. 
This form has an A/D or input quantizer, Q^, digital-to-analog (D/A) or 
output quantizer Q 2 and one internal feedback quantizer Q^. The transfer 
functions from each quantizer to the output are 


T x (z) = D(z) 

T 2 (z) = 1 

b.z -1 + b,z -2 

t 3 (z) " - — , " T o * 

1 + b^z + b 2 z L 

The integrands for equation (4-18) are thus 

V*) T^l/z) (a 0 z 2 + a x z + a 2 ) (a Q + a^ + a 2 z 2 )/b 2 
2 z(z 2 + b^z + b 2 ) (z 2 + b^z/b 2 + l/b 2 ) 

T 2 (z) T 2 (1/z) 


(6-2) 


(6-3) 
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Fig 26. The Direct Form. 
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T 3 (z) T 3 (1/z) (bjZ + b 2 )(b 1 z + b 2 z 2 )/b 2 

z z(z 2 + b 3 z + b 2 )(z 2 + b 3 z/b 2 + l/b 2 ) 

For programming and impulse testing the difference equations for the 
direct form are 


e ± (k) q " e i< k ) + n l (k > 

e Q ( k ) • a 0 e i^^q^ + a l e i^ k “ *)q 3 + a 2 e i^ k “ 2 ^q^ 

“bje 0 (k - Dq^ - b 2 e Q (k - l)q 3 

e o (k) q 2 " e o (k) + n 2 (k) 
e o (k) q 3 " e o (k) + n 3 (k>> 


The filter output variable is e Q (k) q2 . This completes the description of 
the direct programming form. 

For all eleven programming forms the standard notation of for 
the filter input quantizer and Q 2 for the filter output quantizer has been 
adopted for convenience. The transfer functions T^(z) and T 2 (z) are then 
always to be for the input and output quantizers respectively. These 
transfer functions will be identical for all the programming forms as 
given in equation (.6-2) . 

Modified Direct Form 

The modified direct programming form for equation (6-1) is shown in 
Fig. 27. This form differs from the direct form only in the feedback 



Fig. 27. The Modified Direct Form. 
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loop. This form has two internal quantizers; Q 3 is identical to the 
direct form hence T-j(z) is given by equation (6-2) ; 0^ has been added and 
its transfer function to the filter output is displayed below: 


V z > 


1 + b-L 



+ b2^ 


-1 


(6-5) 


The integrand for equation (4-18) for Q 4 becomes 

T 4 (z) T 4 (1/z) z 2 /b 2 

z z(z^ + b^z + b 2 ) (z 2 + b^z/b^ + l/b 2 ) 

For programming and impulse testing the difference equations for 
the modified direct programming form are 


e i (k) q = e i (k) + n^k) 

e 0 (k) = a oei (k) q + a 1 e.(k - l) q + a 2ei (k - 2) q 
+ m(k - l) q 

e o (k) q 2 = 6 ° (k) + n 2 (k) (6 " 7) 

e o (k) q 3 = e o (k) + n 3 (k) 

m(k) = -b ieo (k) q3 - b 2 e o (k - l) q3 

m(k) q = m(k) + n^(k) . 
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Standard Form 

The standard programming form for equation (6-1) is presented in Fig. 
28, This form has two internal quantizers, and Q^. Their transfer 
functions to the filter output are 


V z > - T-r 1 


z + b^z + b 2 


z + b-, 

V 2 > - -i 1 — 

z z + b^z + \>2 


( 6 - 8 ) 


The integrands for equation (4-18) are 


T 3 (z) T 3 (1/z) 


z^/b. 


z(z^ + b^z + b 2 )(z^ + b 1 z/b 0 + l/b 2 ) 


r ,u 2 


(6-9) 


T 4 (z) T 4 (1/z) (z + b 1 )(z + b lZ 2 )/b 2 


z(z^ + b^z + b 2 ) (z^ + b 3 z/b 2 + l/b 3 ) 


The difference equations for this form are 


e ± (k) q = e i (k) + n-^k) 
e o (k) = a O e i (k) q + m 2 ( k " 1 >q 


e 0 ( k ) q " e Q (k) + n 2^ k ) 


m 1 (k) = ct 2 e ± (k)q ~ b^m^k - l) q - b 2 m 2 (k - 1) ( 


( 6 - 10 ) 
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m lOOq = m 1 (k) + n^OO 
m 2 (k) = o^e^k^ + n^Ck - l) q 
m 2 (k)q = m 2 (k) + n^(k) 

where 

“1 = a l " a O b l 

a 2 ~ a 2 “ a 0 b 2 ” b l a l* 

Modified Standard Form 

Again the modified standard form is for D(z) as expressed in equation 
(6-7) and is demonstrated in Fig. 29. This programming form differs from 
the standard form in its feedback loops. The same internal quantizers 
are present as before with a fifth quantizer added. The transfer func- 
tions for the three quantizers are 

T 3 (z) = T 3 (z) in (6-2) 

T^(z) = T 3 (z) in (6-8) 

T 5 (z) = T 4 (z) in (6-5) 

Hence, the integrands for equation (4-18) have been previously shown. 
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For programming, etc. , the difference equations for the modified 
standard form are 

e i (k) q = e^Ck) + n^k) 
e o (k) = aQe^k^ + m 2 (k - l) q 

e o (k) q 2 “ e o (k) + n 2 (k) 

e 0 (k)q 3 " e Q (k) + n 3 (k) 

m^k) = a 2 e i (k) q - b 2 e Q (k) q ^ 

m l( k >q = m i^ k ) + n 4^ k ) 
m 2 ( k ) = m i( k — l)q 

m 2 (k) q = m 2 (k) + n^(k). 

Canonical Form 

The block diagram for the canonical programming form limited to 
the second-order module of equation (6-1) is shown in Fig. 30. This 
form has only one quantizer Q 3 whose transfer function to the filter 
output is given by 


( 6 - 11 ) 


b l e o <k> q. 


T 3 (z) = D(z) . 




nical Form. 
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Therefore, Q 3 has the same effect as the input quantizer on the filter 
output. The difference equations including quantization are shown 
below: 

e i^ k )q = e i^ + n l^ k ) 

m(k) = e i (k) q - b^k - l) q - b 2 m(k - 2) q 
m(k)_ = m(k) + n~(k) 

q • (6-12) 
e D (k) = a Q m(k) q + a^k - l) q + a 2 m(k - 2) q 

e Q (k) q = e Q (k) + n 2 (k). 


Modified Canonical Form 

The modified canonical programming form for the second-order D(z) 
of equation (6-1) is depicted in the block diagram of Fig. 31. This pro- 
gramming form differs from the canonical form by its forward transfer 
paths. By moving the multiplier coefficient from m(k) q to e i (k) q the 
transfer function for is changed: 


T 3 (z) 


a-^z + a 2 


z L + b-^z + \>2 


where 


(6-13) 


“l ' a l ' a o b i 

a 2 “ a 2 “ a 0 b 2 # 
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The integrand for equation (4-18) for this transfer function is 

2 

T 3 (z) T 3 (1/z) (c^z + a 2 )(a 1 z + a 2 z ) 

— 5 ~ • (6—14) 

z z(z + b^z + b 2 ) (z^ + b^z/b 2 + l/b 2 ) 

The difference equations for the modified canonical programming form are 
shown below: 

efOOq = e i (k) + n^k) 

e 0 (k) = a Q e i (k) q + a 1 m(k - l) q + a 2 m(k - 2) q 

e G ( k )q e Q (k) + n 2 (k) (6-15) 

m(k) = e.(k) q - b x m(k - l) q - b 2 m(k - 2) q 
m(k) q = m(k) + n 3 (k). 

The six programming forms discussed to this point have all required 
the programming coefficients and b^ of equation (6-1) , or were easily 
calculated from them. The last five forms which are to be presented now 
will require more effort to find the correct form for D(z) and the proper 
programming parameters for the difference equations. 
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Parallel Form 

The general block diagram for the parallel programming form for 
a second-order D(z) is shown in Fig. 32. The form may be used if 
and only if the second-order module has real poles p-^ and p 2 . Hence, 

D(z) must have the form 

D(z) = a 0 + + ^2 (6-16) 

z-Pl z _ P 2 • 


The constants and R 2 are real numbers representing the residues of 
poles p^ and P 2 , and p^ should be different from P 2 » 

The coefficients g A shown in Fig. 32 must satisfy the following 
relationships : 

glg 2 - R i 

8 3 8 4 = R 2 


(6-17) 


In order to minimize the magnitude of the parameters g^ the following 
choices were arbitrarily made: 

8 i = 

g 2 = R i/gi 

(6-18) 

g 3 = ^I r 2 I 

g 4 — r 2 / ®3 * 


The transfer functions from the internal quantizers to the filter 
output were obtained: 


T 3 (z) 


8 2 

z-p x 


(6-19) 
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Fig. 32* The Parallel Form. 
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g 4 

T (z) = 

z-P 2 • 

The integrands for equation (4-18) corresponding to (6-19) are 

T 3 (z)T 3 (1/z) = -gj z/p 3 

z z(z-pi) (z-l/p,) 

1 ( 6 - 20 ) 

T 4 (z)T 4 (1/z) = -g 4 z/p 2 

z z(z-p 2 ) (z-l/p 2 ) . 

The difference equations for the parallel fora are 
e i(k)q = e i( k ) + n l(k) 

e 0 (k) = a 0 e i (k) q + g 2 m 1 (k-l) q + g 4 m 2 (k-l) q 
e 0 (k) q = e 0 (k) + n 2 (k) 

( 6 - 21 ) 

m l (k) = 8 1 e i ( k ) q + P^Ck-Dq 
m 1 (k) q = m 1 (k) + n 3 (k) 
m 2 00 = g 3 e i(k) q + p 2 m 2 (k-l) q 
m 2 (k)q = m 2 (k) + n 4 (k) . 

Please note that the parallel fora can realize only real poles, but 
it is capable of realizing either real or complex zeroes. 

Cascade Form 

The cascade programming fora for a second-order digital filter 
module essentially factors the module into first order stages and 
realizes each stage individually. If each first order stage is 
implemented in the manner of Fig. 30; the resulting cascade form is 
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shown in Fig. 33. A requirement for this form is that 

D(z) = a Q (z-qj) (z-q 2 ) 

(z-p-^) (z-p 2 ) » 

where q. and p are real zeroes and poles. Also, the following 
i i 

relationships must be satisfied! 


a o = h 8 2% 

g 3 =-g x g 2 (6-23) 

85 "*"""^2^4 

The cascade form has two internal quantizers which are described 
by the transfer functions 
T,(z) = D(z) / g. 

3 1 (6-24) 


T 4 (z) = g 4 z-q 2 
z-p 2 


The integrands for equation (4-^18) are 
T (z) T (1/z) = 1 D(z) D(l/ 2 ) 

_3 17 

Z (6-25) 

T^(z) T 4 (1 /z) = -g 4 (z-q 2 ) Cl-q 2 z)/p 2 
z z(z-p 2 ) (z-l/p 2 ) 

The difference equations for this cascade form are displayed below: 
e (k) = e. (k) + n. (k) 

i 9 1 1 (6-26) 

m x (k) = g 1 e i (k) q + P^Ck-Dq 


m^j^CkJq = m^Ck) + n 3 (k) 
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Fig. 33. The Cascade Form. 
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n» 2 (k) = g^Ck)^ + gg m ^ Ck-l)^ + p 2 m 2 (k-l)^ 
m 2 (k)^ = m 2 (k) + n^(k) 

e 0 (k) = g 4 m 2 (k) q + S5 m 2^ k_1) q 
e Q (k)q = e 0 (k) + n 2 (k) . 

The parameters g is i=l, 5 in equation (6-26) must be found using (6-23). 
Since there are three equations with five unknowns, an arbitrary choice 
for g^ and g 2 is made as follows: 

gl = i.o 

1 (6-27 

82 = /RJ. 

If a Q is zero, this form cannot be realized. 

This completes the cascade form. In summary, this programming form 
is applicable to a second-order digital filter module when it is 
possible to cascade first-order stages programmed in the canonical form. 


Modified Cascade Form 

Of the many possible ways of implementing first-order stages, 
one other technique was selected which employs the modified canonical 
form for each first-order section (see Fig. 34). This programming 
form is labeled modified cascade; it requires D(z) to be factorable 
into real poles and zeroes as in equation (6-22). 

The transfer functions from the three internal quantizers to the 
filter output are given below: 


T 3 (z) = g 3 g 4 (z-q 2 ) 
(z-pq) ( z — p 2 ) 


(6-28) 
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Fig. 34. The Modified Cascade Form. 
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z-q 2 

V Z) " g 4 

z-p 2 

T 5 (Z) = 8 6 

Z- P 2 > 

where the parameters are restricted by 
g l g 4 = a o 

S 2 g 3 = ( Pl _q l )g l (6 " 29) 

g 5 g 6 = ^ p 2 _q 2^ S 4 * 

Since there are three equations and six unknowns, arbitrary choices 
are again made for g^, g 2 , and g^ as follows: 

g i = 

g 2 = (Pj+D/2 

g 3 = (P i -qi)g 1 /g2 (6-30) 

g 4 = a o/gi 

g 5 = ( p 2 +1 )/ 2 

gg = (P 2 -q 2 )g4/g5 

Using these parameters, the following difference equations may be used 
to implement this programming form: 


e i (k) q = e^(k) + n^Ck) 
m 2 (k) = g 1 e i ( k ) q + g 3 mi(k-l) q 
m2(k)q = m 2 (k) + n 4 (k) 

e o 00 = g 4 ® 2 00 + g6 m 3 Ck-1) q 
e 0 (k) q = e 0 (k) + n 2 (k) 
n^OO = g 2 e i (k) q + p 1 m 1 (k-l) q 
m^k^ = m 3 (k) + n 3 (k) 
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m 3 (k) = S 5 m 2^q + P2 m 3^ k-1) q 
m 3 (k) q = m 3 (k) + n^(k) . 

XI Structure 

The last two programming forms to be presented are designed 

for a second-order D(z) with complex poles. The appropriate 

expression for the transfer function is 

D(z) - a 0 + + A? (6-32) 

z+p z+p , 

where a 0 has been previously defined, p and p are complex conjugate 
poles, and A and A are complex conjugate residues. 

The first implementation of (6-32) is depicted in the block diagram 
of Fig. 35. The parameters indicated in the figure are defined below: 

g x = -Re (p) 
g 2 = Im (p) 

(6-33) 

g 3 = 2 Im (A) 
g 4 = 2 Re (A) . 

The two internal quantizers, and Q^, are described by the 
transfer functions 


T,(z) = §2 (6-34) 

zH-b^z+b^ 


V z) = z ~ §1 

zH-b^z+b^ 


The corresponding integrands for equation (4-18) are 


T 3 (z) T 3 (1/z) = g^z 2 /b 2 

z z(z2+b 3 z+b2) (z^+b-^z/b2+l/b2) (6-35) 
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T, (z) T^Cl/z) = (z-g 1 ) (z-g 1 z 2 )/b 2 


z(z 2 +b^z+b 2 ) (z^+b^z/b 2 +l/b 2 ) 


The difference equations for the XI structure are enumerated below: 
e ± (k) q = e i (k) + n^k) 
e 0 (k) = a 0 e ± (k) q + m 2 (k-l) q 
e 0 (k) q = e 0 (k) + n 2 (k) 

m-^k) - g 3 e i (k) q + g^Ck-l) - g 2 m 2 (k-l) q 
m l ^k> q = m l00 + n 3 (k) 

n^GO = g 4 e i^ k) q + Sl m 2^ k-1 )q + ^l^^q 
m 2 (k) q = m 2 GO + n^(k) . 


X2 Structure 

The last programming form presented in this paper is the X2 
structure of Fig. 36. The transfer function D(z) must be expressed 

in the format of equation (6-32) in order to use this form. 

This programming form has two internal quantizers whose transfer 
functions to the filter output are 


t 3 (z) = §3 Z + (§ 2 8 4 " Si g 3> 


+ b^z + b2 


T 4 (z) = g^z - (g 2 g 3 + gjg 4 ) 
o 

z + b^z + b 2 

where 

g 1 =-Re (p) 
g 2 = Im (p) 
g 3 =~Im (A) 
g 4 = Re (A) 


(6-37) 


(6-38) 
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The integrands for equation (4-18) are 

T 3 (z)T 3 (1/z) - (g 3 z+6 x ) (g 3 z+6 1 z 2 )/b 2 

z z(z 2 +b-jz+b 2 ) (z 2 +b 3 z/b 2 + l/b 2 ) 

T 4 (z)T 4 ( 1 /z) = (g 4 z+6 2 ) (g 4 z+6 2 z 2 )/b 2 

z z (z^+b-^z+bp (z 2 +b^z/b 2 +l/b 2 ) 

where 

6 1 = 8 2 8 4 " S l 8 3 
<$2 =_ §2 s 3 “ g l s 4 • 


The difference equations for the X2 structure are listed below: 
e i (k) q = e^k) + ri 1 (k) 

e 0 (k) = a 0 e i (k) q + g 3 m 1 (k-l)q + g 4 m 2 (k-l)q 
e 0 (k) q = e 0 (k) + n 2 (k) 

®2. (k) = g-jD^ (k-l) q - g 2 m 2 (k-l) q 
m l(k)q = m lOO + n 3 (k) 

m 2 (k) = 2 ei(k) q + g 1 m 2 (k-l) q + g 2 m 1 (k-l) q 

m 2^ k ^q “ m 2 (k) + n 4 (k) . 

This completes the X2 structure. 

Summary of Programming Forms 

This section has summarized the essential characteristics of 
eleven programming forms for a second-order digital filter module. 
All of the equations necessary to perform steady-state, statistical, 
and error bound analyses have been determined. A pattern may be 
observed in the formats of the relations for equation (4-18) , the 
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residue evaluation equation for statistical analysis. All of the 
integrands fall into the following formats: 


F-^z) = Y3(y 0 2 2 +y 1 2+Y 2 ^ (Y 0 +Yl z ~hr2Z 2 ) 

z(z 2 +b^z+b 2 ) (z 2 +biz/b 2 +l/b 2 ) 


(6-41) 


or 

F 2 (z) = Y2Cz-Y^) (1-Yiz) 

z(z-Y 0 ) (z-1/y 0 ) • 

Table 2 displays the respective equations for each programming form 
using equation (6-41) . 

Many other characteristics of each programming form should also be 
investigated; for example, the coefficient sensitivity and the deadband 
effects are also important for good digital filter operation. 
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Table 2. Integrands for (4-18). 



Modified 

Standard 

q 3 

Q 5 

F 1 

F 1 

F 1 

0 

0 

0 

j 

b l 

0 

0 

b 2 

1 

1 

l/b 2 

l/b 2 

l/b 2 

Canonical 

q 3 

F 1 

a o 

a l 

a 2 

l/b 2 

Modified 

Q 3 

F 1 

0 

1 

al 

«2 

i/ b 2 

Canonical 

1 

! 






Parallel 

Q 3 

F 2 

Pi 

0 

-*2 /p l 



Q 4 

F 2 

p 2 

i 

0 

2/ 1 
-S4 /p l 

1 

Cascade 

^3 

F 1 

a o 

a l 

”1 

a 2 

1/g l b 2 


Q 4 

f 2 

1 

p 2 

^2 

4 /p 2 

I 

! 

Modified 

Q 3 

F 1 

■ 

0 

1 


" 1 
g 3 g 4 /p l p 2 

Cascade 








^4 

F 2 

P 2 

q 2 

-84^2 



^5 

f 2 

P 2 

0 

- 4 '* 2 
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Table 2. Integrands for (4-18). (Cont'd) 


Programming 

Form 

Quantizer 

Format 

Parameters 

Y o 

a 

^2 

^3 

XI Structure 

q 3 

F 

1 

a 

0 

1 

8 2 /b 2 


^4 

F 1 

0 

1 

_8 1 

i/b 2 

X2 Structure 

Q 3 

F 1 

0 

g 3 

<5 X 

l/b 2 


Q4 

F 1 

0 

g 4 

6 2 

l/b 2 























VII. COMPUTER AIDED DESIGN 


In the design of complex system, the digital computer serves as 
an essential tool in synthesis and design verification. Computer 
aided design (CAD) programs are effectively employed in the synthesis 
of digital filters in three ways: transfer function synthesis, 

coefficient quantization, and programming form selection. 

Transfer Function Synthesis 

The digital computer has been used extensively in the design of 
digital filter transfer functions [30,71-74]. Nonrecursive designs 
using linear programming has been implemented by Rabiner [73] while 
Parks and McClellan [72] using polynomial interpolation techniques. 
Rabiner et al [30] also used a steepest descent technique to obtain 
FIR filters with minimax error in selected bands. 

Recursive digital filters have been synthesized using sampled 
data transformations by Golden [71]. Robinson and Robinson [74] 
have demonstrated a CAD program for taking z-transforms. Steiglitz 
[75] has used nonlinear optimization techniques to obtain recursive 
digital filter approximations to arbitrary frequency responses. 
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Coefficient Quantization 

Avenhaus [48] has investigated the effects of coefficient optimi- 
zation for reducing the coefficient wordlength. A given filter is 
designed and its coefficients are founded. Then an optimizing search is 
undertaken to find other sets of coefficients which meet the design 
criteria with a shorter wordlength. 

Much work is left to be done in the proper quantization of digital 
filter coefficients and CAD will surely play a major role in future 
developments in this area. 

Programming Form Selection 

A CAD program, listed in [49], has been developed which analyzes 
the signal amplitude quantization errors in the eleven programming 
forms presented in Chapter 6. The program, written in FORTRAN IV, 
is an aid to implementing digital filters for any application, the only 
restriction is that the filters be expressable as second order stages 
as shown in equation (6-1). 

General 

The filter implementation program actually consists of eleven 
parts, one for each programming form discussed in the previous section. 
Each programming form is analyzed using steady-state, statistical, 
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and upper bound techniques. The system weighting constants K . , K ct . • 

SSJ SlJ 

and K ubi are calculate d using the equations of Table 1. K and K 

ssj stj 

were computed by both equations for debugging purposes; K ub j was 
determined using the second equation. A weighted average of these 
constants was also used: 

^aj = H^ssj + ^2^stj + ^3^ubj * (7-1) 


where 


^1 + ^2 + ^3 = 


Therefore, a weighted average error can be calculated by 


[el = l K .h! 
o wa ^ waj 3 


(7-2) 


The weighting constants may be adjusted by the designer to emphasize 
the steady state, statistical, or upper bound errors. 

The filter implementation program may be used in two modes, one 
for stored-program computers and one for special-purpose hardware; the 
two modes are distinguished by the manner in which the quantizer step 
lengths are chosen. In both the assumption is made that truncation 
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(or LSB-1) quantizers are used in the system. All errors must be 
halved by the user if roundoff quantizers are present. 

I 

In the stored program mode the maximum error hj of the jth 
quantizer is fixed by first simulating the ideal digital filter 
response to a "worst case" step input, which is an A/D input word 
of all "ones." During the transient response to this step, the 
maximum value of the filter output and internal variables is recorded. 
After the simulation has run a sufficient number of iterations for conver- 
gence, sav 100, the maximum values are rounded up to the nearest 
power of two. Since the computer wordlength is a fixed number L r , 
the quantizer intervals are found by 



(7-3) 


h^ is always assumed equal one. 

In the special-purpose computer mode the register lengths are 

f 

not fixed; therefore, a different method is used to find hj. The 
philosophy of this mode is to balance the effect of each quantizer 
in the system so that they all have relatively equal error contribu- 
tions. This balancing is done by dividing equation (7- 2 > *val 
(with h| = 1) : 


CeJ 


wa 


= 1+1 


Kwal 



(7-4) 


Each term in the summation is forced to be less than or equal one 
(to insure that the A/D will introduce an error as large or larger 
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than the other quantizers) by choosing 

f 

hi 1 Kwal * (7-5) 

Kwaj 

f 

A further restriction is that hj be a power of two; hence the ratio 
Kwaf/Kwaj is rounded down to the nearest power of two to find the 
actual hj to be used; 

hj = [Kygi/Kyaj] rounded down* (7-6) 

Flow Charts 

Fig. 37 demonstrates the flow of information in the main section 
of the filter implementation program. The input data to be given to 
the program is summarized below: 

1) Transfer function coefficients in (.6-1) 

a o » a i* a z» ^1* ^2 

2) Register lengths 

A/D, d/A, and wordlength of the stored-program computer 
or coefficient wordlength for the special-purpose computer. 

3) Weighting coefficients in (7-1) 

This is all the information needed to completely analyze the quantiza- 
tion errors for all the programming forms. 

The first major task in the program is to find the poles and 
zeroes of the transfer function D(z) and to set three flags which omit 
those programming forms which are unrealizable. The main program 
then calls a subprogram for each realizable programming form. Each 
called subprogram completely analyzes the quantization errors 
characteristic to that particular form and prints their detailed 
description. At the end of the program, final summaries of each 
program mode are listed for easy cross-reference. 
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Fig. 37. Flow Chart of Main Program 
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A general flow chart describing a subroutine for any given 
programming form is shown in Fig. 38. The first task is to calculate 
all of the parameters needed for the difference equations of the 
specified programming form; next, these parameters are quantized. 

The simulation difference equations are then calculated once for 
the step response and once for each quantizer in the system. During 
these simulations the system constants Kg S j , K s tj , and K^j are 
calculated. Finally the steady-state, statistical, and upper bound 
errors are calculated, as well as the weighted average error of 
equation (7-2), for both modes of program operation. 


Source Listinc 


The filter implementation program consists of approximately 1800 
source statements and is available in £49]. Also, a limited 
number of printed listings are available from Auburn University. 


Summary 


The filter implementation program has been developed using an 
IBM 360/50 using FORTRAN IV and OS360. In its final form the program 
takes approximately 3.5 minutes to compile and 25 seconds to load and 


execute. The execution time may be trimmed by limiting the simulation 
iterations to a smaller number, say 10 to 20. 

Now the CAD program will be used to analyze two digital filters, 
one for each program operating mode. 
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Fig. 38. Flow Chart for Programming Form 
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Stored Program Mode 

Consider the second order digital filter 

D ( Z ) = £ ± • 75z ..' f -0.-il 5 . (7-7) 

z 2 + ,50z + .0525 

Suppose that this filter is to be realized using a 16-bit minicomputer 
using a 11-bit A/D and 13-bit D/A as input-output equipment. The com- 
puter-aided design (CAD) program may be used in the stored-program mode 
of operation to aid the designer in programming the minicomputer. Table 
3 is the final summary of quantization errors attributed to the filter 
above for its realizable programming forms. The D(z) in (7-7) has real 

poles and zeroes; therefore, the XI and X2 structures may not be used. 

Using the weighted average errors in Table 3, the CAD program rec- 
ommends that the filter in (7-7) be programmed by first the modified canon- 
ical form; and second, the parallel form. Note that all the programming 
forms give relatively good results; this is due to the fact that the 
internal quantizers and output quantizers contribute only a minor part 
of the total quantizing error. The A/D and D/A wordlengths chosen in 
this example are responsible for these results. 



TABLE 3. Final Summary: Stored-Program Computer Mode 
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Special-Purpose Computer Mode 

Suppose that a special-purpose computer is to be constructed to im- 
plement the following z -domain transfer function: 

D( 2 ) - z 2 ~ 1.862z + .895 # (7-8) 

z 2 - .2500 

Again, if an 11-bit A/D is to be used, the CAD program gives the results 
shown in Table 4. From the table, the weighted average error suggests 
that the direct form is best; the modified canonical form, second. 

Direct form . The program prints out an analysis of each programming 
form which may be used for (7-8). See Table 5. The system error weight- 
ing constants (K gg , K gt , and K^) are summarized as well as the X’s of 
(7-1), the maximum quantizing error h' of each quantizer, and the form 
factor. The form factor is interpreted as follows 

FORM = I,J, 


where 


I « total number of bits for the register 
J = number of bits to the right of the binary point. 

A negative J indicates the least significant bit has a value greater 
than one. From Table 5, h^ = 2h^ and h^ * 4h^. The CAD program always 
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TABLE 5: The Direct 


STEADY-STATE ANALYSIS 
KSS(l) = 0.04 A 
KSS(2) » 1.000 
KSS(3) » 0.333 

STATISTICAL ANALYSIS 
KST(l) - 1.426 
KST(2) = 0.577 
KST(3) = 0.149 

ERROR BOUND ANALYSIS 
KUB(l) = 5.009 
KUB(2) » 1.000 
KUB(3) * 0.333 

SPECIAL-PURPOSE COMPUTER MODE 

LAMBDA(l) » 0.333 H(l) - 1.0 

LAMBDA (2) = 0.333 H(2) - 2.0 

LAMBDA (3) ■ 0.333 H(3) - 4.0 

STEADY-STATE ERROR - 3.377 
PERCENT Q1 = 1.3 
PERCENT Q2 - 59.2 
PERCENT Q3 - 39.5 

RMS ERROR = 3.177 

PERCENT Q1 = 44.9 
PERCENT Q2 = 36.3 
PERCENT Q3 = 18.8 

MAXIMUM ERROR BOUND » 8.343 
PERCENT Q1 - 60.0 
PERCENT Q2 * 24.0 
PERCENT Q3 - 16.0 

WEIGHTED AVERAGE ERROR = 4.966 
PERCENT Q1 = 43.5 
PERCENT Q2 - 34.6 
PERCENT Q3 = 21.9 


Printout 


FORM - 11,0 
FORM - 10,-1 
FORM » 9,-2 
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assumes - 1. Also, the form factor of Q 2 , the output quantizer, sug- 
gests that a 10-bit D/A may be used. 

Modified canonical form . The CAD program output for the modified 

f j I 

canonical form is shown in Table 6. Note that h^ » 2h^ and hj * h^ for 
this programming form. 

Closed-Loop Comparison 

The second-order digital filter in (7 _ 8) has been analyzed in [53J 
for a closed-loop sampled-data control system. The block diagram for 
the control loop is shown in Fig. 39. Statistical and upper bound tech- 
niques were employed to design the compensator of the control loop for 
both the direct and modified canonical forms; system simulations were 
employed to verify the results. Table 7 presents a comparison of the 
open-loop results of this paper and the closed-loop results of [53]* 

Note that they agree very closely. 

One observation should be made at this point. The register lengths 
determined by the open-loop design procedures of this paper are in gen- 
eral larger than those required in closed-loop applications. Stable 
feedback systems generally tend to reduce the maximum values of the dig- 
ital filter variables and thus the number of bits needed to represent 
these variables in the special-purpose computer. 
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TABLE 6: The Modified Canonical Printout 


STEADY-STATE ANALYSIS 
KSS(l) - 0.044 
KSS(2) = 1.000 
KSS(3) - -0.956 

STATISTICAL ANALYSIS 
KST(l) - 1.426 
KST(2) - 0.577 
KST(3) « 1.303 

ERROR BOUND ANALYSIS 
KUB(l) - 5.009 
KUB(2) » 1.000 
KUB(3) - 4.009 

SPECIAL-PURPOSE COMPUTER MODE 

LAMBDA (1) - 0.333 H(l) - 1.0 FORM - 11,0 

LAMBDA(2) * 0.333 H(2) - 2.0 FORM - 10,-1 

LAMBDA(3) - 0.333 H(3) - 1.0 FORM « 12,0 

STEADY-STATE ERROR - 3.000 
PERCENT Q1 » 1.4 
PERCENT Q2 = 66.7 
PERCENT Q3 - 31.9 

RMS ERROR = 3.884 

PERCENT Q1 = 36.7 
PERCENT Q2 * 29.7 
PERCENT Q3 = 33.6 

MAXIMUM ERROR BOUND « 11.019 
PERCENT Q1 » 45.5 
PERCENT Q2 - 18.1 
PERCENT Q3 - 36.4 

WEIGHTED AVERAGE ERROR - 5.968 
PERCENT Q1 ■ 36.1 
PERCENT Q2 = 28.8 
PERCENT Q3 *= 35.1 




Fig* 39. Control-loop [53]. 
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TABLE 7: 

Open-Loop Versus 

Closed-Loop 

Programming 

Form 

Open-Loop 

Results 

Closed-Loop 
Results [17] 

Direct 

h£ - 2hJ 

h’ 


hj - 4h^ 

h 3 - 4h ; 

Modified 

h 2 ' 2h i 


Canonical 

h 3- h i 

h^ - .5hi 
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Conclusion 

This section has presented a computer-aided design technique useful 

in implementing digital filters expressed as z-domain transfer functions. 
Two examples have been given to illustrate the stored-program and special- 
purpose modes of operation of the CAD program. Also, the prograro.which 
analyzes the filter's "open-loop" quantization errors, gives results closely 
matching a "closed-loop" design. This CAD program should be used as a 
tool for obtaining a "first guess" at the best way to program a digital 
filter. If a closed-loop simulation is available for the system in which 
the digital filter will be used, then the CAD program design may be ad- 
justed to give better loop performance. 

Although the program as presented has been designed for second-order 
modules, it can be used as a subroutine in larger programs to match pole- 
zero pairs for higher order realizations, or to indicate the proper 
cascade ordering of second-order modules. The CAD program may be a power- 
ful tool to the digital filter (or controller) designer if its results 
are properly interpreted. 



VIII. APPLICATIONS OF DIGITAL FILTERING 


Digital Filtering has found many diverse applications in recent 
years. This section lists several of them and points the interested 
reader to the open literature for detailed descriptions. 

The following list presents typical applications for digital filters 

1. Sampled-Data Control Systems 

a. General [38, 77] 

b. Pendulous Integrating Gyroscopic Accelerometer [78, 79]. 

c. Saturn V Thrust Vector Control [80, 81] 

2. Speech Processing 

a. General [82] 

b. Vocoder [83] 

c. Equalizers [84] 

3. Radar and Sonar Signal Processing 

a. General [85, 86] 

b. MTI Filters [87, 88] 

c. Tracking Filters [89, 90] 

4. Spectral Analysis and Synthesis 

a. Narrow Band Filters 

b. FFT [91] 

c. Frequency Synthesis [92] 

5. Vibrations and Acoustic Testing [93] 

6. Image Processing 

a. General [24, 94, 95] 

b. Image Enhancement [94, 95, 96] 

c. Pattern Recognition [97] 

7. Seismic Processing [7, 9] 

8. Biomedical Processing [94, 95, 97, 98] 

9. Synthesis of Speech and Music [99] 

Many other applications of digital filtering are also important 
with the number of new ones ever increasing. 
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I . INTRODUCTION 


In recent years a trend has been developing to replace analog sys- 
tems with digital systems. This rate of replacement has been directly 
related to the technological advances in the manufacturing of digital 
logic building blocks. With the advent of large-scale- integration, a 
particular class of digital networks, called digital filters, has be- 
come economically practical in such areas as stabilization of control 
systems, spectrum analysis, voice and' speech analysis, radar, medical 
electronics and virtually any other analog filter function 11,2]. 

Digital filtering has been defined in PARTS ONE and TWO as a 
computational process consisting of digital multiplications, additions 
and delays whereby one sequence of numbers is transformed into another 
sequence. This transformation may be specified by a transfer function 
in the z-domain, D(z) , or by a set of linear difference equations 
with constant coefficients. Assuming knowledge of these coefficients, 
digital filter realization procedures [1,3, 4, 5] consist of the design 
of a digital system to solve these difference equations. The difference 
equations may be solved with a software program and a general purpose 
computer or with the use of a special-purpose (SP) computer [6, 7, 8, 9, 
10,11,12], a technique which has become increasingly popular. 
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In Che SP computer realizations, a particular digital filter pro- 
gramming form is selected and the computer is designed accordingly 
{13,14,15,16,17], Particular attention must be given to assure that 
the hardware organization meets the system specifications for coefficient 
quantization, signal amplitude quantization, and quantization noise 
levels introduced into the system by th.e digital filter implementation. 

At the present time no systematic design procedure has been developed 
to accomplish these goals. Typically one designer specifies the digital 
filter coefficients and another specifies a hardware implementation. 

The state-of-the-art in digital filter implementation is represented 
in {1,3,7,8,9,10,12,18,19,20,21]. Perhaps the most interesting are the 
IC model in [1] and the programmable design of 112]. The IC model is 
available from Autonetics Division of North American Rockwell. A 
digital filter implemented with this technology is small and can 
realize third-order filters at sampling rates of up. to 5kHz. However, 
poles must be real and the parallel programming fora is the only one 
available. The commercial units also have restricted programming 
forms, or implementation is done by frequency transformations which 
limit their use to applications in which minimum time delay and high 
speed sampling are not specified. It has been shown in [14] that some 
of the programming forms have different characteristics, and it is 
desirable in many cases to be able to select the programming fora. 

The need for a selectable programming form along with the desirable 
features of LSI implementation offer a challenge to the system designer. 
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Add the necessity for real-time fault diagnosis and standardized CAD 
procedures to the list and the system design goals are complete. 

Now that the theory of digital filtering has been presented in 
PART TWO, we will examine four mechanization techniques for digital 
filters. The four techniques are 1) general-purpose computers, 

2) mini-computers, 3) spqcial-purpose computers, and 4) FFT hardware. 
A discussion of all techniuqes will be presented starting with 
mechanization (implementation) by general-purpose (GP) computers. 




II. GENERAL PURPOSE COMPUTER IMPLEMENTATIONS 


Of the four implementation techniques, the GP computer is the 
least attractive, with the main reason being that most GP computers 
possess excessive computing capabilities to be used only for differ- 
ence equation calculations. If this were done, there would be large 
portions of the computer hardware that would never be used thereby 
making this type implementation overly expensive. 

GP computers do have a useful application in digital filter 
implementation in that they may be used to simulate other implementation 
designs (an example being by special-purpose computers) or for real- 
time programming of a GP computer to implement a digital filter as well 
as other computational chores. Let us now look at these two aspects 
of using a GP computer in the design of a digital filter system. 

Simulation 

The most common implementation of a digital filter is by special- 
purpose computer. When designing a special-purpose computer for the 
implementation of the filter in a particular programming form, one of 
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the first steps that must be done is deciding on word length require- 
ments for the input word, output word and internal variable (m(kT - T), 
m(kT - 2T), etc. of the difference eqs.) wordlengths and possibly 
arithmetic schemes. This can be accomplished by techniques such as 
the CAD program presented earlier. Once a design is recommended, it 
is good engineering practice to simulate the system on a GP computer 
to verify all the design parameters. With most higher level languages, 
logical programming may be done such that every aspect of the design 
may be simulated. If this approach is taken the system designer may 
"change something" and observe its effects; this technique may be 
used to "optimize" the final system design. 

As an example of a digital filter implementation simulation, the 
program below was written in FORTRAN and run on an IBM 360 digital 
computer to simulate the "range-switching" filter described in [8] 
employed in a nulling type control loop. The program was written so 
that the "range-switching" effects on the output of the loop could be 
observed and the effect the wordlengths had on the output response for 
a particular input, which in this case is a sine wave of specified 
amplitude and frequency. 
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FORTRAN SOURCE PROGRAM FOR SIMULATION OF PIGA LOOP 
WITH A NOISE INPUT 


SOURCE DECK 

C TIME DOMAIN SIMULATION OF THE COMPENSATED SYSTEM 
DIMENSION X1(2),X2(2),X3(2), D(2) 

COMMON/COM2/RM(1001) ,ITER,TES 
COMMON/ C0M1/XP(3,2) ,BQ(1) ,C(2) 

1 FORMAT (1H ,13»2X,1P9E13.4) 

CALL INPUT 

2 FORMAT (1H1) 

N=1 

H=l. 0E+05 
WN=184. 0 
T=0. 001 
GP=56. 2 
GT^321000.0 

C GDA IS THE TOTAL LOOP GAIN 

GDA=0. 083 
W1=SIN(WN*T)/WN 
W2=COS (WN*T) 

W3= (1. 0-W2) /WN**2 
W4=(T-W1)/H 
W5= (1 . 0-W2) /H 
W6=WN*SIN (WN*T) 

W7=W6/H 

GDIG=GAD/GT/GP/ (T-Wl) *H 
WRITE (6 ,10) GDIG 
10 FORMAT (1H0,7HKDIG = .1PE11.4) 

TEST=+2000. 0*980. 0/4. 0 
DO 669 NAGL-1,2 

669 XP (NAGL, 1)=0. 0 
8 BDC=0. 0 
BQ(1)=0. 0 
X1(1)=0.0 
X2 (1)=0. 0 
X3 (1)=0. 0 
D(1)=0. 0 
C(1)=0.0 
COFS=0. 0 
DO 5 1=1, ITER 
R=TES*980. *RM(I) 

4 WRITE (6,1) I,BDC,BQ(N) ,X1(N) ,XP(1,1) ,D(N) ,R,COFS ,C(1) 

C BEGIN ANALOG PORTION SIMULATION 

XI (N+1)=X1 (N)+W1*X2 (N)+W3*X3(N)+W4*(R-GR*D(N) ) 

X2 (N+1)=W2*X2 (N)+W1*X39N)+W5*(R-GR*D(N) ) 

X3 (N+l ) =W2 *X3 (N) +W6 *X2 (N) +W7 * (R-ST*D (N) ) 

BDC=GP*X1 (N+l) 
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END ANALOG PORTION SIMULATION 
BEGIN DIGITAL UNIT SIMULATION 
UI=BDC 

CALL DIGCOM (UI,YP) 
D(N+1)=GDIG*YP 
C END DIGITAL UNIT SIMULATION 

COFS=GT*D(N+l) 

XI (N)=X1 (N+l) 

X2 (N)=X2 (N+l) 

X3 (N)=X3 (N+l) 

5 D(N)=D(N+1) 

STOP 

END 

SUBROUTINE DIGCOM(Ul ,YP) 
COMMON/COM1 /SP (3 , 2 ) ,BQ(1) ,C(2) 
A0=1.0 
Al=-119./64. 

A2=57./64. 

B1=0.0 

B2=0.0 

FX=256. 

UI=UI/3.0*FX 
EX=1. 0 

CALL ROUND (UI,EX,FX) 

UP-UI 

Z BQ(1)=UP*3. 0/FX 

UI=UI*3.0/FX 
IF(ABS(UP)-16.0) 2,3,3 

2 C(2)=0. 0 
GO TO 4 

3 C(2)=1.0 
UP=UP/16. 

IUP=UP 

UP=IUP 

FX=16. 

4 IF(C(2)-C(1) ) 4,6,7 

5 XP(1,1)=16. *XP(1,1) 

XP(2,1)=16. *XP(2,1) 

AX=4. 0 

BX=63. 75 

CALL ROUND (XP (1, 1) , AX, BX) 

CALL ROUND(XP(2 ,1) ,AX,BX) 

GO TO 6 

7 XP(1,1)-XP(1,1)/16.0 

XP(2,1)=SP(2,1)/16. 

AX=4. 0 
BX=63. 75 

CALL ROUND (XP(1,1) ,AX,BX) 

CALL ROUND (XP(2,1) ,AX,BX) 

6 SP (1 ,2)=-Bl*XP (1,1)-B2*XP (2 ,1)+UP 
XP(2,2)=SP(1,1) 
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AX=4. 0 
BX=63. 75 

CALL R0UND(XP(1,2) ,AX,BX) 
YP=(A1-AO*B1)*XP(1,1)+(A2-AO*B2)*XP(2,1) 
1 +AO*UP 

CX=64. 

DX=255./64. 

CALL ROUND (YP,CX,DX) 

DO 1 1=1,2 
1 XP(I,1)=XP(I,2) 

C(1)=C(2) 

YP=YP/FX*3. 0 

RETURN 

END 

SUBROUTINE INPUT 
C 

C RANDOM INPUT 

COMMON/ COM2/RM(1001) ,ITER,TES 
TES=200. 

ITER=301 

NRANB=6 

CALL RANBIT(NRANB) 

CALL RCON1 (35187269) 

RMAX= (2 . **NRANB-1 . ) /2 . 

D019 1=1, ITER 
RM(K)=IRAN(5) 

19 RM(I) = (RM(I)-RMAX) /RMAX 
RETURN 
END 

SUBROUTINE ROUND (A, AN, BN) 

X=ABS (A) 

S=A/X 
IX=X*AN 
XQ=IX 
XQ=XA/AN 
IF(XQ-BN) 1,2,2 

1 A=S*XQ 
RETURN 

2 A=S*BN 
RETURN 
END 
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Real Time Programming 

A digital filter implemented on a general-purpose computer, 
whether large or small, is said to be realized by real-time pro- 
gramming. The machine language version (translated from some higher 
level language) of the difference equations must execute quickly 
enough to meet the sampling rates imposed by the system specifications 
In some applications the general-purpose computer will handle other 
calculations as well and will be "time-shared" to perform both duties. 
Other times a small process control computer can be dedicated solely 
to the digital filter calculations. An example system is shown in 
Fig. 2.1. 

Generally speaking, future trends will be to design special- 
purpose computers to shoulder the digital signal processing tasks, and 
relieve the general-purpose computer for more complicated tasks which 
exploit its entire computational power as embodied by its versatile 
instruction set. 
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Fig, 2,1. A GP computer being used as a digital 
filter in a discrete control loop. 



Ill . MINICOMPUTER IMPLEMENTATIONS 


A minicomputer “Implementation of a digital filter as described 
in [22] will be discussed. Only one reference is used as a background 
since it is the only one that has been seen in the literature of digital 
filtering. It will be sufficient since any other minicomputer imple- 
mentation would follow the guidelines presented. 

Hardware Requirements . 

The hardware used for the minicomputer implementation is shown in 
Fig. 3.1. It consists of a Honeywell H316 minicomputer with two 4096- 
word memory modules, a 10— bit analog— to— digital (A/D) converter, a 12— bit 
digital- to-analog (D/A) converter, a crystal-controlled real-time clock 
and the ASR-33 teletypewriter. 

H316 minicomputer . The H316 is a GP minicomputer with a 16-bit 
wordlength. Arithmetic is performed in two registers, A and B, and 
it is a one-address machine with the A register serving as the accumulator 
which will be described in detail later. The memory is divided into 
sectors or pages of 512 words each, with the computer having the 
capability to reference any of the 512 words within a certain sector. 
Single-level indexing and/or multiple-level indirect addressing can be 
used to address words outside the current sector or the base sector. 

With respect to the arithmetic instructions of the computers 
instruction set, there are two modes of operation: single precision 

and double precision. Each mode of operation may be entered by the use 
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Fig. 3.1. Hardware used in minicomputer implementation. 
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of one instruction, "SGL" for single-precision arithmetic operations 
and "DBL" for double precision arithmetic operations. When operating 
in the single-precision mode, the A register is used solely as the 
accumulator. It is 16-bits long with the left-most bit being the 
sign bit and the 15-bits to the right being the most-significant 
through the least significant of the magnitude bits which are in a 
two’s complement code. When operating in the double precision mode, 
the A and B registers are used as the accumulator with the sign bit 
being in the left-most bit position of the A register. The rest of 
the A register contains the 15 most significant bits of the double 
precision word with the 15 least significant bits being contained in 
the 15 right most bit positions of the 16 bit B register. The left 
most bit position of the B register does not take part in arithemtic 
operations. 

When performing the ’’add” instruction which will have to be done 
many times in difference equations calculations, the contents of the 
addressed memory word are added to the contents of A leaving the sum 
in A for single-precision addition. If done in double-precision the 
contents of the addressed memory word (two memory locations for double 
precision) are added to the contents of the A and B registers and the 
sum left in them. 

The same procedure occurs for multiply for the single or double 
precision mode. The addressed word in memory is multiplied by the word 
stored in the A or A and B registers and the product left in the A or 


A and B registers. 
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Fixed point arithemtic is used for all difference equation calcu- 
lations. Since an imaginary binary point is assumed, after a multipli- 
cation instruction is executed, the computer shifts the product as 
required to align the binary point. 

Input to the minicomputer is accomplished through 16 input bus 
lines into the A register. Several peripheral devices may be connected 
to this bus as inputs to the computer. In the case of the, minicomputer 
implementation of a digital filter this bus inputs information from the 
A/D, ASR-33, and the real-time clock. 

Output is accomplished through 16 output bus lines which are tied 
directly to the A register and always reflect its contents. For the 
digital filter implementation, the output device is the D/A or the ASR-33. 

The different input devices are checked by the computer by placing 
a code unique to each device on the address bus. 

Most peripheral devices are slower than the computer, thereby 
making the computer spend much of its time waiting for a peripheral 
device to perform its function. It is for this reason that it is 
practical to let the computer process other information while a particu- 
lar peripheral is performing its I/O function. Then when the peripheral 
is finished, it can inform the computer and the computer can give it 
another command. 

The method of informing the computer of the completion of a task 
is called an interrupt. When a peripheral interrupts the H316, it 



3-15 


finishes executing the instruction presently being performed and then 
performs a subroutine jump indirectly through a dedicated memory 
location. In short, the dedicated memory location contains the address 
of a subroutine to which the computer jumps when an interrupt occurs. 
Within this subroutine the computer may poll the peripherals to find 
out which one interrupted. 

An A/D converter is used as the input interface element to the 
computer. The A/D which was interfaced to the H316 is a bipolar con- 
verter having a range of -lOv to +10v. It has a 10-bit plus sign-bit 
output which is input into the most significant 10-bits of the A 
register. 

The D/A converter accepts and transforms the binary output of the 
computer into an analog voltage. A hold register is employed so that 
the output voltage will remain constant until the next output occurs. 

The D/A used in the minicomputer implementation was built from Honeywell 
y-Pac DTL logic. The converter is built from three cascaded Honeywell 
CE-071 four-bit converters which consist of a resistive ladder plus 
switching network. 

To provide for a sampling rate other than that determined by the 
computers execution time, a real-time clock was employed. It initiates 
each cycle consisting of input, calculation, and output and is built 
as a peripheral which furnishes periodic (sample rate) interrupts to 
the computer. When an interrupt occurs, the computer goes through one 
cycle and then waits for the next interrupt before it goes through the 
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cycle again. This allows the operator of the minicomputer to obtain 
any desired sample rate. 

Operating System 

It is the purpose of the minicomputer implementation to be able 
to realize in real-time one of eleven different digital filter programming 
forms. It is the function of the operating system to set up, control, 
and possibly run diagnostic tests if something goes wrong, on the 
minicomputer and its peripherals. 

A functional block diagram of the operating system (OS) is shown 
in Fig. 3.2. Solid arrows indicate a passing of control from one 
routine to another, while dotted arrows indicate a passing of parameters. 
Only one filter form is shown, but it should be remembered that eleven 
such forms are present with similar links to the operating system 

Briefly, to realize a digital filter, a particular form is picked 
and the parameters which determine the transfer function are input. 

The OS will then type back these parameters if desired. Once the filter 
form is set up, the OS is instructed to begin execution of that form. 

Let us now discuss the different parts of the OS. 

Executive . The executive routine (EXEC) initially types a question 
mark on the teletype. Whenever the question mark appears the operator 
types in one of four commands: MODIFY, LIST, RUN, or TEST. The first 

three refer to a particular programming form and are followed by a 
number between one and eleven. The TEST command refers to one of seven 
diagnositc routines and should be followed by a number from one to seven. 



Diagnositcs 
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Fig* 3.2* Block diagram of operating system. 
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After one of the four commands is typed in EXEC turns control over to 
one of the four routines having the same name. Let us briefly discuss 
these routines. 

1. Modify . The modify routine inputs the coefficients, quantization 
formats, and sample rate for a particular filter form. EXEC determines 
which of the eleven forms has been typed in following the command MODIFY, 
then transfers control to the modify routine, passing the filter form 
number as a parameter. 

2. List . The list routine types out the coefficients of a pro- 
gramming form followed by the quantization formats and finally the 
sample period. 

3. Run . To begin the filter processing the operator would type 
in RUN followed by the number of the form he desires to use. 

RUN has a list of all entry points of the filter forms. When 
the RUN routine is entered it immediately obtains the address of the 
normal entry point and passes it to the interrupt processer which will 
need it at a later time. Then RUN selects the sample period which the 
user has specified for that filter form and outputs it to the RTC. 

Next, RUN sets the mask of the real-time clock (RTC) and teletype, 
starts the RTC, types out a question mark, and transfers control to 
the initialization entry of the filter form specified. The filter form 
makes its first pass and hangs up in the idle routines at the end. 

While in the routine the RTC should interrupt. 



3-19 


Interrupt processor . When the interrupt occurs, control is passed 
to the interrupt processor. This routine must identify what caused 
the interrupt and act accordingly. 

User's interface . The user interface consists of the teletype 
routines plus the data conversion routines. The teletype routines 
are relied upon by all the other routines which have to communicate 
with the user. The teletype routines handling mumerical data rely on 
the conversion routines to convert from decimal to binary and binary 
to decimal. 

Diagnostics . Seven diagnostic routines are implemented in the 
OS to test the hardware and the software structure. One of these routines 
may be executed by typing in the request TEST followed by a number from 
one to seven. The errors that are checked for are divided into three 
categories: hardware errors, errors in the programming of the OS, and 

last, user errors. 

This completes the discussion of the OS. We will now look at the 
assembly programs. 

Assembly Programs . 

Each of the eleven filter programming forms is realized by a 
separate subroutine which has the following format: 

. ENTRY 
INPUT 

CALCULATION 

OUTPUT 

TIME DELAY 

PRECALCUATION 

IDLE 

EXIT 
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There are two entry points to each program: one being an initializa 

tion entry point which zero's the internal variables the first time 
through the program and the second a regular entry point that is 
entered every time except the first. After entering the normal entry 
point a "start A/D" command is given and, while waiting on the input 
to become available, a partial sum is formed. As soon as the input 
arrives it is shifted to a correct format, multiplied by Ag, and the 
sum is then completed. The sum is then quantized for output, presented 
to the D/A, then quantized in a different format for storage and feed- 
back. If overflow is detected during quantization, the word is saturated 
i.e. filled with the largest possible number. 

After the output is complete, the internal variable must be shifted 
to perform time delay. Then the partial sum for the next pass is begun. 
Just enough of the formation of the partial sum is left for the next 
pass to occupy the arithemtic unit while waiting on the A/D. During 
the "idle" period, the RTC interrupts and the interrupt subroutine 
directs control back to the normal entry point. 

The coefficients as well as the three shift instructions used in 
the quantizing routines are declared as external names so that they 
may be altered by the OS. 

As an example of one of the eleven assembly language programs the 
assembly language program for a second order D(z) in modified canonical 
programming form is shown below. 
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SUBR 

MCAN1 ,ENT1 

SUBR 

MCAN2 , ENT2 

ENT 

SHFT61 ,S1 

ENT 

SHFT62 ,S2 

ENT 

SHFT63 ,S3 

ENT 

COEF6.AO 

BEL 


* INITIALIZE INTERNAL VARIABLES 

ENT1 CRA 


STA 

XM1 

STA 

XM2 

★CALCULATE 

OUTPUT DIFFERENCE EQ 

ENT2 OCP 

' 41 

LDA 

XM1 

MPY 

AL1 

DBL 


DST 

TEMP 

LDA 

XM2 

MPY 

AL2 

DAD 

TEMP 

DST 

TEMP 

INA 

'1041 

JMP 

*-l 

SI LRS 

4 

STA 

El 

MPY 

AO 

DAD 

TEMP 

SGL 


STA 

SGN 

S2 LLS 

9 

SSC 


JMP 

OKI 

LDA 

SGN 

CSA 


LSA 

='77777 

SRC 


TCA 


OKI OTA 

'40 

JMP 

*-l 


CALCULATE FEEDBACK DIFFERENCE EQ. 
LDA El 
MPY ONE 
DBL 

DST TEMP 
LDA XM1 
TCA 
MPY 


START A/D 


INPUT FROM A/D 
WAIT FOR INPUT 


B1 
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DAD 

TEMP 


DST 

TEMP 


LDA 

XM2 


TCA 



MPY 

B2 


DAD 

TEMP 


SGL 



STA 

SGN 

S3 

LLS 

3 


SSC 



JMP 

OK2 


LDA 

SGN 


CSA 



LDA 

='77777 


SRC 



TCA 


* PERFORM 

TIME DELAY 

0K2 

STA 

XM 


LDA 

XM1 


STA 

SM2 


LDA 

XM 


STA 

XM1 


ENB 



NOP 



JMP 

*-l 

XM1 

DBP 

0 

EO 

BSS 

1 

XM 

BSS 

1 

El 

BSS 

2 

XM2 

BSS 

2 

TEMP 

BSS 

2 

SGN 

BSS 

1 

AO 

OCT 

10000 

AL1 

OCT 

-22753 

AL2 

OCT 

7357 

B1 

OCT 

3146 

B2 

OCT 

231 

ONE 

OCT 

10000 


END 



Experimental Results 

Experimental results were obtained of the minicomputer implementation 


previously described 
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First it realized the transfer function of a Euler integrator 


D(z) 



(II-l) 


in the direct form at its maximum sampling rate (5.5 KHz). The response 
was obtained for an input square wave and as wished, the output was a 
triangular waveform with a fine-grained stair-stepped appearance. 

Secondly it realized the transfer function of a digital differentiator 

D(z) =1 - z -1 (II-2) 


in the direct programming form. It's response to a triangular wave- 
form, a square wave, was as expected. 

Lastly, the transfer function of a digital oscillator was realized 

D(z) = - (II-3) 

1 - 2cos(2irfT)z _1 + z" 2 


where T is the reciprocal of the sample rate (5.5 KHz) and f is the 
frequency of oscillation. The equation was programmed with b-^ =-1.75 
which resulted in an output frequency of 450 Hz as predicted by 
Eq. (II-3) . 

Experimentation demonstrated that the direct and canonical forms 
had the highest maximum sampling rate. The direct, canonical and 
modified canonical forms have minimum sample intervals of less than 200 
y-secs. The modified direct form has a minimum interval of about 225 y-secs. 
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The parallel and cascade forms have a minimum interval of about 275 p-secs, 
while the remainder of the forms have intervals of approximately 
300 p-secs. 

The number of instructions required to implement the filter forms 
(including coefficient storage), ranges from 98 for the canonical to 
109 for the modified cascade. The entire operating system occupies 
approximately 5000 memory locations including indirect links in the 
base sector. 

This concludes the discussion on minicomputer implementations of 
digital filters. From this discussion it was seen that a minicomputer 
can be adapted well for a real-time digital filter implementation; in 
fact, much better than the larger SP computers because of the smaller 
size. We will now look at even a smaller digital computer implementation, 
that of implementation by special-purpose computers. 



XY. SPECIAL-PURPOSE COMPUTERS 


It is obvious to one that for most realizations of a digital filter, 
the general purpose computer and the minicomputer approach have several 
disadvantages. The most easily seen disadvantage, as previously mentioned, 
is the wasted hardware incurred because of the relative simplicity of 
the difference equations that must be calculated for a realization. It 
is for this reason that the special-purpose computer approach to realiza- 
tion is taken for a majority of the applications of digital filtering. 

It will be shown in the following discussion of special-purpose (SP) 
computer realizations that they are the most economical (hardware wise) 

and demonstrate a great amount of versatility. 

The realization techniques by SP computers will begin with the very 

earliest method, which was sample and hold devices with analog networks 
and conclude with present day commercial models that are available on 
the market. 

Implementation by Sample and Hold Devices with Analog Networks . 

The first SP computer realizations of a digital filter were by sample 
and hold devices with analog networks \2 3] . There are two main ways of 
realizing a digital filter in this manner with them being: 1) series 

discrete data networks, and 2) feedback discrete data networks. There 
is a third way of realization which is a combination of the above two 
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that will not be discussed since it is felt it is not essential to 
illustrate the realization technique. 

The series discrete data network which is used for the realization 
of the discrete transfer function of a digital filter is shown in 
Figure 4.1. The transfer functions of the two systems in Fig. 4.1 
are related by 


D(z) " Vc (!) 


(4-1) 


_Ts 

Since the transfer function of the zero-order hold is (1-e )/s,.it 

can be obtained from Eq. (4-1) that 


21 . 1 ■ —^T D(z > (4-2) 

1 - z 

From this it is seen that given a specific D(z), the transfer function 
G c (s) of the discrete-data network can be determined from Eq. (III-2) 
by taking the inverse z-transf orm. 

If G c (s) is to be an RC realizable transfer function, all the poles 
of G c (s) must be simple and lie on the negative real axis of the s-plane 
with the exception of the origin and infinity. The zeroes of G c (s) may 
be located anywhere in the s-plane. Therefore, G c (s)/s can be expanded 
into the following form by partial fraction expansion: 


G (s) 
c 

s 



\ 

s + s k 


(4-3) 




(b) Equivalent series discrete-data network. 


Fig. 4.1 








3-28 


where A q and ^ are constants and s k = [k = 1, 2 , 3, * * * , m ] are simple 
negative real poles. The z-transform of Eq. (tu-3 ) is 


z[ 


G c (s) 



+ 1 - 

k=l 

1 





(4-4) 


which has simple positive real poles inside the unit circle |z| =1, 
with only one pole at * - 1. Comparing Eq. (4-4) with Eq. (4-2), 
it is seen that in order for G c ( s) to represent an RC network, the 
discrete transfer function D(!z) must have the following properties: 

(1) The number of poles of D(z) must be equal to or greater than 
the number of zeroes of D(z). 

(2) The zeroes of D(z) are arbitrary in location. 

(3) The poles of D(z) must be simple, real and positive, and 
lie inside the unit circle |z| = 1 in the z-plane. 

It can be shown that for a feedback discrete-data network, the 
feedback structure with a zero-order hold and the RC network shown in 
Fig. 4.26 is equivalent to the digital filter of Fig. 4.2a [23]. 
transfer function of the two systems are related by 


D(z) 


1 

1 + G ho H(z) 


z [ Slsl 1 1 r 1 - D(z) . 

1-z-l 1 DC*) 1 


or 


(4-5) 
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(a) Digital filter. 



RC Network 


(b) s Equivalent feedback discrete-data network. 


Fig. 4.2 
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In order to realize H(s) by an RC network, the discrete transfer 
function D(z) must have the following properties: 

(1) D(z) must have the same number of poles and zeroes. 

(2) The poles of D(z) are arbitrary. 

(3) The zeroes of D(z) must be simple, real, positive, and lie 
inside the unit circle of the z-plane. 

The reasons behind these restrictions are discussed in detail in 
123] and therefore are not repeated here. 

In summary it may be said that the listed restrictions and limited 
flexibility of the above implementation techniques make their use and 
practicality almost negligible. 

Hybrid Implementation . 

After observing the difference equations of several of the pro- 
gramming forms which may be used for the realization of a digital 
filter, such as the direct and canonical forms, one can see that these 
equations might be computed by a device which performs the arithmetic 
functions of addition, subtraction, multiplication and time delay. 

This suggests the use of summing amplifiers, constant multipliers and 
delay elements for the implementation. Also, knowing that digital 
filters may be implemented by SP and GP computers suggests the realization 
of a digital filter by hybrid techniques, with hybrid meaning that analog 
and digital methods are used for the implementation 110] . For most hybrid 
realizations digital techniques are used for analog-to-digital (A/D) 
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conversion, digital-to-analog (D/A) conversion and time delay with analog 
techniques used for the arithmetic functions of addition, subtraction, 
and multiplication. This type realization exploits the best and most 
natural functions of both analog and digital elements to eliminate a 
majority of the restrictions placed on the realization of the previous 
section of the literature which used sample and hold devices with analog 
networks. Another advantage of the hybrid realization which will be . 
illustrated shortly is that once a D(z) is obtained, the filter can be 
realized directly from it. This will eliminate the need for additional 
mathematical manipulation required for the derivation of an s-plane 
transfer function from the D(z) as was required by the previous implemen- 
tation technique . 

The hybrid implementation of a digital filter fits itself to a 
majority of the programming forms of a given D(z). For simplicity sake, 
we will look at the hybrid realization of a D(z) in two programming 
forms; the direct form and the canonical form. These two forms usually 
have the simpliest difference equations which must be implemented for 
their realization. 

For most hybrid designs an A/D converter functions as the sampling 
device as well as an interfacing element for the input to the filter. 
Integrated circuit buffer registers are usually used to store previous 
values (in digital form) of an intermediate variable that is internal to 
the controller. D/A converters provide the interface for the output of 
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the filter. Variable resistors at the input of a summing amplifier are 
usually used to adjust the coefficients of the compensation function 
continuously over a wide range of values. The equivalent of a zero-order 
hold device is realized at the output of the filter as a result of the 
digital data-storage elements within the hybrid unit. 

The hybrid filter is designed, in this case, to realize almost 
any compensation function up to three zero's over three poles in the 
z domain with any sampling frequency up to several thousand hertz. If 
a filter of order higher than three is required, several second or third 
order filters are cascaded until the desired order is obtained. This 
is usually done over constructing a higher order filter because of the 

greater coefficient sensitivity for a higher order filter. It is seen that 
the hybrid controller is extremely versatile and can be used to realize 
a wide range of sampled— data compensation functions; thus, as previously 
mentioned, it is not necessary to redesign the unit to change the 
compensation function. 

To obtain an insight into the design of a hybrid digital filter 
for a given D(z), a hybrid implementation will be derived for the direct 
and canonical programing forms. 

Let us look at the direct form realization first. A second order 
transfer function for a realizable digital filter can be written as 


D(z) = 


, -1 t -2 

a o + a l Z + a 2 Z 

-1 -2 
1 + b^z + b£Z 


E (z) 
E t (z) 


» 


( 4 - 6 ) 
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where the coefficients a 1 and are real numbers and z = e . 

Equation (4-6) can be rewritten as 

E q (z) = a E (z) + a^z ^(z) + a 2 z 2 E 1 <z) - b^^z ^(z) - 

b 2 z" 2 E (z) C 4 ’ 7 ) 

4 O 

or in the time domain as the difference equation 

e Q (kT) = a Q e 1 (kT) + aje (kT - T) + a^OcT - 2T) - b^CkT - T) - 

b 2 e Q (kT - 2T) C4-8) 

where f = 1/T is the sampling frequency. Fig. 4,3 is a block diagram 
s 

of the hybrid realization of the second order direct programming form 
with double lines denoting digital information and the single lines 
analog information. A detailed explanation of the components used for 
the realization will be given after the block diagram of the canonical 
form is given. 

In order to realize the canonical form, we have previously seen 
that for a second order D(z), the following difference equations must 
be realized. 

m(kT) = e^kT) - b^kX - T) - b 2 m(kT - 2T) (4-9) 

e (kT) = a n m(kX) + a,m(kT - T) + a.m(kX - 2T) 

O U 1 4 


(4-10) 
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The block diagram that results for a canonical hybrid realization is 
shown in Fig. 4.4. It is easily seen ho w this form results after 
careful consideration of the form of the two difference equations that 
must be implemented. 

The actual construction of the canonical implementation will now 
be discussed in more detail. This form was chosen to be discussed 
because of several advantages it has over the direct programming form. 

One of the primary reasons, which is obvious from Figs. 4 .3 and 4.4, 
is that the canonical form is much more economical. For example, the 
direct form requires two A/D converters whereas the canonical requires 
one. Also, if n is the order of the numerator of the D(z) being realized 
and m the denominator, the direct form requires n + m delays whereas the 
canonical requires the greater of n and m. This is also the case for 
D/A converters. It can be said that in general anytime a D(z) is to 
be realized by a hybrid realization, choose the programming form which 
requires the least hardware. 

Let us now discuss the analog and digital components of a hybrid 
implementation. Observing the canonical realization of Fig. 4.4, the 
first step in a hybrid realization is the conversion of analog informa- 
tion into digital information by the use of an A/D. There are two types 
of A/D's which might be used, the first being the successive approximation 
type and the second the count-up- to type. 

Fig. 4.5a is a block diagram of a successive approximation type 
A/D converter. The comparator compares the A/D input signal m Q (the 

Preceding page blank 
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LSB -*■ Least significant bit 
MSB -*• Most significant bit 


(a) Successive approximation type 
A/D converter. 



(b) Counter type A/D converter. 


Fig. 4.5. A/D converters. 
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signal m Q is the output of the feedback summing amplifier as shown in 

Fig. 4.4) to the quantized output m 0( j of a D/A converter. The signal 
m oq is determined wholly by the digital word m(kT) which is in sign 
magnitude code. The logic circuitry is programmed to "search for" or 
"home-in on" the analog signal m 0 . Successive approximation type 
converters are faster (commercial models are available that will convert 
an analog signal to an 8 bit sign magnitude code approximation in 
200 nsec) than the counter types which will be described shortly, 
therefore they are used in a majority of applications. 

Fig. 4.5b is a block diagram of the counter type A/D converter 
In general it is slower than the successive approximation type converter 
and is therefore used when the input analog signal m Q is of lower 
frequencies. This type of converter operates on the principle of 
letting a binary counter count until its output decoded through a D/A 
converter is equal to the input signal. When this occurs the counter 
ceases operation and its output bit sequence (m(kT)) at this time is 
a sign magnitude code approximation of the analog input m 0 . Before the 
converter can be ready for the next conversion the binary counter will 
have to be set such that all its bits are logic "0", with this being 
completed before each conversion. 

Of the two types of converters discussed above, it is recommended 
that the successive approximation type be used for digital filter 
implementations because of its ability to handle high frequency input 
signals. 
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The time delay indicated by Eqs. (4-9) and (4-10) can be 
provided by clocked flip-flops. This means the information at the 
inputs of the flip-flop is held or stored until the clock terminal 
is pulsed, at which time the information that is on the input is 
transferred to the output terminals. 

The D/A converters shown in Fig. 4.4 are the ladder type networks 
commonly used in constructing D/A f s, From the Figure it is seen 
that the contents of the time delays constitute previous values of 
m(kT) . These digital words are decoded by the D/A converters into 
analog signals and are then available for further analog processing; 
that is, multiplication and summation. 

From Eqs. C4-9) and (4-10) it is seen that the a^ and b^ are 
real numbers and may be positive or negative. This suggests the use of 
operational amplifier circuits to perform the arithmetic operations of 
multiplication and summation. From Fig. 4.6, which illustrates in 
detail the summation and multiplication techniques for a second order 
implementation in canonical programming form, it is seen that the 
algebraic sign of the coefficients is obtained by inverting the output 
m^ of the i^ D/A (multiplying by -1) so that m-^q and -m.^ are available. 
From Fig. III-6 it is also seen that the magnitude of the coefficients 
is obtained by variable input resistors to an operation amplifier. 

Any hybrid digital filter implementation should generally be 
implemented in the same procedure as the above. If there are variations, 
they are usually small, and are left up to the individual designer. 
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Digital Implementation . 

The adaptation of a small SP computer for the implementation of 
digital filters only seems natural after a careful consideration of the 
requirements that must be met for the realization of the difference 
operations of a particular programming form. Let us consider the 
requirements of the difference equations of an arbitrary programming 
form since the requirements would hold for all forms. Let our choice 
be the difference equations required for the realization of a D(z) 
in the modified canonical programming form. These equations for a 
second order D(z) are shown below: 

e 0 (kT) = a 0 e ± (kT) + ct^mCkT - T) + <x 2 m(kT - 2T) (4-11) 

m(kT) - e^kT) - b^kT - T) - b 2 m(kT - 2T) (4-12) 

Observing the above equations we see that for them to be physically 
realized a device must be used which can add, subtract, multiply, 
perform time delay, truncate and provide data storage. A device which 
can do all of this is a small SP computer. All of the components of a 
digital computer can be organized into four main functioned units as 
shown in Fig. 4.7. 

Considering the functional requirements of a digital filter, we see 
that the Arithmetic Unit of the computer can perform the addition, 
subtraction, multiplication, and the truncation required for the 
realization of the digital filter. The Memory can be used to perform 
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Fig. 4.7. Four functional units of a digital computer. 
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the time delay, coefficient storage, internal variable storage, and 
input/output storage required for a filter implementation. The Input/ 
Output functional block of the digital computer will accomodate the 
A/D and D/A converters if required for interface for a digital filter 
realization. The Control Unit functional block of the digital computer 
will accomodate the controller for the digital filter and the data 
transfer logic which insures correct routing of data for the real- 
ization of the required difference equations. At this point we see 
that all of the arithmetic, storage, input/output, and control require- 
ments necessary for the implementation of a digital filter have all 
been incorporated into the four functional blocks of a digital computer; 
i.e., a digital filter may be realized by a small SP computer and its 
functional diagram is shown in Fig. 4.8 [7]. The computer is said 
to be small because it will be designed to only calculate the difference 
equations for a particular programming form. Another reason the resulting 
SP computer is small is because of reduced word lengths that are required 
thereby requiring less hardware. 

Now that it has been shown that a small SP digital computer can 
be used for the realization of a digital filter, it is now appropriate 
to discuss how one would go about designing a SP computer realization 
of a digital filter and the considerations that must be made while doing 
this. 

The design consists of three parts: first, the determination of 

quantization levels in the computer (input quantizing, round-off errors, 
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Fig. 4.8. Functional diagram of a digital filter. 
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and filter coefficient quantizing); second, the logical design of the 
computers components; and third, the interconnection of the computers 
components to implement the above mentioned quantization levels . 

The first step in the design procedure, the determination of 
quantization levels, will not be covered here since it would be a 
reiteration of earlier sections of this work. As a reminder 
though, one could use several techniques to do this, among them being 
the CAD program and different quantization error analysis techniques. 

Next consider the second step, which has not been presented 

previously and will be discussed in depth. 

The second step in the design of any SP computer implementation of 

a filter, as previously mentioned, consists of the logical design of 
the computer components. All of the components will be grouped into 
four main component groups with these being, 1) Input/Output components, 
2) Arithmetic components, 3) Memory components and 4) Controller 
components. The discussion will begin with the design of the input/ 
output components. 

1. Input/Output Components 

When designing the input/output equipment for a digital filter the 
first decision to be made is that of what type information will be the 
input and output of the filter, i.e. is the input/output of the filter 
going to be in analog or digital form. 

If the input and output of the filter is in digital form the design 
problem will usually be minimal since the filter normally operates on 
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digital inputs and outputs in digital form. The only problem that might 
be encountered if the application of the filter is such that it will 
have digital input/output is that of synchronization between the device 
that is supplying the filter its input and the filter. Provisions must 
be made when designing the filter such that it will accept a digital 
input word from a device which is supplying it. 

For many applications of digital filtering, the 
filter will be operating in an environment composed mainly of analog 
signals which will necessitate input/output interface elements for the 
filter. These interface elements will be A/D converters for the input 
to the filter and D/A converters for the output of the filter as shown 
in Fig. 4.9. 

First let us discuss the selection of an A/D converter. The first 
decision to be made is how fast should its conversion speed be. In 
general this is dictated by the frequency of the input signal or the maxi- 
mum sample rate of the filter. If the input analog signal is of high 
frequencies and the filter is to sample at a high rate the successive 
approximation type converter previously discussed usually would be 
better since it is faster than the counter type. Next, a 
consideration of the number of bits required for the digital approxi- 
mation of the analog signal would be required. This can be determined 

by knowing the maximum analog input voltage (V ) to the filter and 

max 

the maximum allowable quantization step length h. If these two parameters 
are known the number of quantization steps required can be determined by 
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Fig. 4.9. A digital filter with its interface elements. 
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V 

Number of Quantization Steps = max = QS (4-13) 

h 

Knowing QS, the number of magnitude bits (n) required for the A/D 
converter can be chosen if n is the smallest integer such that 

(2 n - 1) >. QS (4-14) 

Now that the bit length and speed of the A/D converter is known, 
the next step is the construction or selection of a commercial A/D. 

In selecting a commercial A/D, there are many distributors available. 

All one has to do to find them is to thumb through an electronics 
or computer oriented magazine. A few pointers to remember when 
selecting commercial A/D converters are that the faster their conversion 
speed, the more they cost, and if internal reference voltages are supplied 
they are also more costly than when the user supplies the references. 

If a user were to construct his own A/D converter, its cost would 
also be determined by its bit length, speed, and the construction of 
circuits to supply reference voltages if they are not already available. 
There is abundant material available on how to construct successive 
approximation and counter type converters. 

Chosing the output interface element, the D/A converter, is generally 
one of the easier design tasks of designing a digital filter. If a D/A 
converter is required, to construct an analog approximation 
of the digital output of a filter, there are two basic types from 
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which to choose. One type has a current output proportional to the magni- 
tude of the digital word it is converting and the other type has a voltage 
output. Which type is chosen depends on the application of the digital 
filter being designed. Most D/A converters are constructed from resistive 
latter type networks. 

There are many commercial types of D/A's on the market today and 
their manufacturers are easily found by simply, as for A/D's, thumbing 
through computer oriented magazines. 

There are several characteristics of D/A's that must be considered 
when buying one. One consideration is if the D/A will have to be unipolar 
or bipolar. Some D/A's will only operate in the unipolar or bipolar mode, 
and others can be connected to operate in both. Another important 
characteristic to consider is the speed in which the conversion is made 
and the settling time of the D/A. If operation of the digital filter is 
in the low frequency range, not as much attention will have to be paid 
to this as if it were operating in a high frequency range. As anyone 
might speculate, the faster the conversion time and settling rate, the 
higher the price. 

2. Arithmetic Unit 

When designing the arithmetic unit of a digital filter, the first 
major decision that has to be made is that of parallel arithmetic operations 
or serial arithmetic operations. Both schemes have been proposed and both 
have their advantages and disadvantages. The question as to which method 
is better can only be determined by the designer and his use for his filter. 



3-50 


In general parallel arithmetic operations are used when there is a 
desire for speed and serial arithmetic operations are desirable for a 
minimum of hardware. Techniques of performing parallel and serial arith- 
metic operations will be considered in more detail shortly. 

The second decision that must be made when designing an arithmetic 
unit of a digital filter is that of what type binary code to use. Since 
the arithmetic unit must be able to add, subtract, multiply and perform 
truncation, it would seem logical to use a signed one's complement or a 
signed two's complement binary code. The reason for using these codes 
over a straight signed magnitude code is because, with them, subtraction 
can be performed by an addition process. Also, multiplication can be 
performed by a shifting and addition process using these codes. Of the 
signed one's complement code and the signed two's complement code it 
seems that a majority of the time the two's complement code is used. 

Now an example of an arithmetic unit organized in a parallel and 
serial fashion will be presented, starting with a parallel organization. 

This particular arithmetic unit is able to add, subtract, multiply, 
and perform truncation. The code used by the computer is a signed 
two's complement code. 

The arithmetic unit performs the addition and subtraction process 
of two n bit words, A and B, by use of n full-binary-adders (FBA) 
connected in the configuration shown in Fig. 4.10. If addition is to 
be performed (A + B), A and B are applied to the input of the adder circuit 
and the output will be the sum of A and B. If the arithmetic operation 
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Fig. 4.10. Full-binary-adders used for addition and 
subtraction of two's complement binary 
numbers. 




3-52 


(A - B) is performed, then A is applied as it is to the input of the 
adder and B is two's complemented and applied to the input. The resulting 
output of the adder is (A - B) . 

Now that it has been demonstrated that we can add and substract 
two words in two's complement code with FBA circuits, an accumulator will 
be defined and it will be shown that with an accumulator and shift 
registers, all the arithmetic operations of addition, subtraction and 
multiplication can be performed. 

A binary accumulator consists of a register, which stores a binary 
number (the augend) in signed-magnitude form and upon receiving another 
binary number (the addend) in the same form adds the second number to the 
first and then stores the sum in the register. The logic diagram of a 
parallel binary accumulator is shown in Fig. 4.11. Each flip-flop 
in the accumulator functions as a modulo 2 counter. The augend is 
initially stored in the accumulator and, during the addition process, 
each flip-flop counts parallel incoming pulses representing the addend 
bits and generates a carry pulse to the next significant bit when the 
flip-flop changes its state from 1 to 0. 

To illustrate the use of the parallel binary accumulator in performing 
the multiplication and summation process , observe the difference equation 


e 0 (kT) = a oejL (kT) + Cl m(kT - T) . 


(4-15) 




Fig. 4.11. Parallel binary accumulator 
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Let 

a 0 = +1.25 = +(1.01) 2 

coefficient 

c^ = +0.50 = +(0.10) 2 + ^ o 

- -*■ 1 

e^kT) = +6. = +(110. ) 2 

variable 

m(kT - T) = -3. = -(011.) 2 
then 


a e ± (kT) = +110. c^CkT - T) = -Oil. 

x +1.01 x +0.10 

+110. -000. 

+ 00.0 - 01.1 

+ 1.10 - 0.00 

+ lll.io -001.10 


and e 0 (kT) = + 111.10 - 001.10 = +110.0 = +(6.) 10 . 

A simplified block diagram of a shift register and parallel binary 
accumulator implementation of the example difference equation is included 
in Fig. 4.12. The following sequence of events is offered to explain 
the operation of this implementation. 

1. Set the accumulator output to zero. 

2. Load c^[+0.10] and m(kT - T)[-011.] into the two shift registers. 

3. Apply an accumulate pulse, a shift pulse, another accumulate 
pulse, a second shift pulse, and a third accumulate pulse 
(abbreviate this sequence by asasa; the output of the accumu- 
lator is now c^m(kT - T) [-001.10]. 


4. Clear the shift registers. 
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5. Load aQl+1.01] and e^ (kT) [+110. ] into the two shift registers 
while leaving cjm(kT - T) shared at the accumulator output. 

6. asasa; the accumulator output is now a 0 e-£(kT) 

+ c^CkT - T) [+110. 00]. 

Thus the difference equation is implemented. 

There are other designs for accumulator of Figs. 4.11 and 
4.12. One design which is used frequently uses FBA's and clocked 
J-K flip-flops. This accumulator design is shown in Fig. 4.13. This 
accumulator operates in the same manner as the previously discussed 
design when used in calculating difference equations. The FBA's are 
used to add the new input to the accumulated sum already in the accumulator 
(stored at the output of the JK flip-flops) and once the new sum is 
obtained an accumulate pulse is applied which clocks the J-K flip-flops 
and causes the input bits of the flip-flops to be transferred to the 
outputs where it will be stored until the next accumulate pulse arrives. 

As stated previously, serial arithmetic may also be used for digital 
filter implementation. Let us consider it now. 

Serial arithmetic is used mainly for two reasons: 1) there is 

a savings in hardware which will be demonstrated shortly and 2) serial 
arithmetic provides for an increased modularity and flexibility in the 
digital circuit configurations. 
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Input to accumulator 



Output of accumulator. 


Fig. 4.13. Parallel binary accumulator 
constructed from FBA T s and 
clocked J-K flip-flops. 
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The two's-complement representation of binary numbers is appropriate 
for digital filter implementation using serial arithmetic because addi- 
tions may proceed, starting with the least significant bits, with no 
advance knowledge of the signs or relative magnitudes of the numbers 
being added. 

The serial arithmetic unit, as the parallel implementation, must 
be able to add, subtract, and multiply and truncate. 

The block diagram of a serial adder (subtractor) is shown in 
Fig. 4.14. 

Briefly it operates in the following manner. 

1. All registers and the delay flip-flop are cleared. 

2. Words A and B, which are to be added or subtracted, are 
shifted into the shift registers to the left of the FBA 
with their binary points aligned. The shifting stops when 
the LSB of A or B reaches the LSB position in the register. 

All registers are now properly set to perform the addition 
(subtraction) process. 

3. Now shift all registers to the right and the delay flip-flop 
for the carry at the same time. The proper sum of difference 
will be shifted into the register on the right. 

If subtraction is being performed by the above circuit, the subtra- 
hend will have to be two's complemented before it is loaded into its 
appropriate register. This operation can be performed with a simple 



\ | j l | Shift Register 

——Shift registers 

1 Shift and add (sub) lead 

Addend or subtrahend register 


Fig* 4.14. A serial adder (subtractor) circuit. 
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sequential circuit which, for each input word, passes unchanged all 
initial least significant bits up to and including the first "1" and 
then inverts all succeeding bits. The circuit that will do this is 
depicted in Fig. 4.15 [18] • 

A serial multiplier configuration is shown in Fig. 4.16. This 
multiplier will only multiply two positive two’s complement numbers, 
which does not restrict it, since the circuit of Fig. 4.15 can be 
used to two's complement any negative number that must be used in the 
multiplication process. Also, one's complement numbers may be multiplied 
by this circuit since if they are negative they may be complemented by 
a simple circuit to obtain the positive form. If this is done for a 
two's complement or a one's complement code, the sign bit of the resulting 
product will have to be retained and if it is negative, the resulting 
product term will have to be complemented appropriately. 

Briefly, the serial multiplier of Fig. 4.16 works in the following 
manner. There are three shift registers Z, X, and Y, a serial adder 
and a half adder. The sign bit leads the serial word as indicated by 
the subscripts of the letters Z, X, and Y. All bits except bits X L and 
Y^ of registers X and Y form a combined shift register. Register Z is 
also a circulating register. 

Initially, the multiplier is stored in register Y with the sign bit 
in Yp Next the multiplicand is serially transferred to register Z with 
the sign bit in Z^. Sign bits Z^ and Y-^ do not move during the succeeding 
shifting of the contents of registers. 
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The value of 1 or 0 of the least significant bit of the multiplier 
in Yg determines whether or not the number bits of the multiplicand 
are to be added; and the addition, if there is one, is carried out. 

Also during this addition time, the circulating register Z restores its 
original contents and the partial product is serially inserted into 
register X, occupying bits X£ to Xg. The combined register is then 
shifted 1-bit to the right; during this shift, any carry bit left in 
the delay flip-flop is shifted into X£, and the least significant bit 
is now at Y 2 . The next addition begins at Xg but not at bit Y£. 

After the right shift, the least significant bit of the multiplier 
in Yg is lost , as Y^ is not a part of the combined register. Yg now 
contains the second least significant bit of the multiplier, which 
possibly could initiate another multiplication. 

This process of addition and right shifting continues until all 
multiplier number bits are shifted out of the combined register. After 
this, the product is available in the combined register with the most and 
the least significant halves of the product being stored respectively 
in the X and Y registers. When there is no round-off the sign bit will 
be in Yj. When there is round-off, 1 is inserted into bit Xg, and the 
sign of the product is inserted into bit X-^. After this the number in 
register X is in desired order [24], 

When using the multiplier described above, the product term addition 
required for the completion of the difference equation calculation pro- 
cess can be handled by the serial adder previously described. 
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There are other techniques of parallel and serial multiplication 
that will not be discussed here. One technique which makes parallel 
multiplication much faster is the Wallace technique described in [25]. 
Likewise, there is a serial multiplication technique presented in [18] 
which makes serial multiplication much faster. 

Usually when a difference equation is calculated by a digital 
filter, it is reduced in bit length before it is stored in memory or 
the output register. It is the purpose of the reduction logic of a 
digital filter to do this. Most reduction logic either performs 
truncation or round-off with most being truncation type. 

The circuitry required for reduction logic is usually simple since 
it is usually a combinational logic circuit. Shown in Fig. 3,14 
is reduction logic which will truncate a signed magnitude binary word 
with 14 magnitude bits (8 to the left of the binary point and 6 to the 
right) such that it has 6 bits to the left of the binary point and 2 
to the right. The lower four bits that are truncated are omitted from 
the truncated word and anytime the untruncated word has a weighted bit 
in one of the two most significant bit positions, every bit of the 
truncated word will be a weighted representation, i.e. the truncated 
word is saturated. 

3 . Memory Design . 

The memory of a digital filter is used for coefficient storage, 
interval variable storage which performs time delay, input word storage 


and output word storage. 
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There are two types of memory that are usually used for digital, 
filter implementations: 1) flip-flop type memories, and 2) read-only- 

type memories (ROM) . 

A majority of the time flip-flop type memories are used for internal 
variable storage (m(kT - T) , m(kT - 2T) , • • • ,m(kT - nT)) which is required 
for the implementation of the time delays and also for the storage of 
the input and output variables of the filter. Shown in Fig. 4.18 is 
a flip-flop type memory for the storage of m(kT - T) and m(kT - 2T) 
required for the realization of a second order D(z) in the modified 
canonical programming form* The word lengths in this figure are 8 
magnitude bits with 6 of them to the left of the binary point. For 
this type memory construction a time delay is performed when the J-K 
flip-flops are clocked; m(kT) becomes m(kT - T) and simultaneously 
m(kT - T) becomes m(kT - 2T). The input and output storage registers 
are constructed the same way as the first column of flip-flops in 
Fig. 4.18. 

ROM type memories are used mostly for coefficient storage. The 
coefficient values are loaded only once into the ROM and they are 
read out when they are needed for an arithmetic operation. 

Digital filters may also employ simple single pole, double throw 
switches for coefficient storage. This is usually done for versatility 
as the coefficients may be changed by a simple toggle of a switch. 

The only disadvantage of this is that the coefficients cannot be changed 
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while the filter is operating in real-time. It must be stopped and 
then started over again. The reason for this is that the switches 
can't all be manually thrown during one sample period which would be 
required. 

4. Controller Design . 

The controller of a digital filter ensures the proper calculation 
of the difference equations being realized by the filter. As an 
example of what a controller must do, let us consider the calculation 
of the difference equations of a second - order D(z) in the modified 
canonical programming form. The difference equations are shown below 


e Q (kT) = aQe^(kT) + c^m(kT - T) + c 2 ni(kT - 2T) (4-16) 

m(kT) = ei(kT) - b x m(kT - T) - b 2 m(kT - 2T) (4-17) 

Eq. (4-16) will be calculated first as it is desirable to obtain 
an output as soon as the input is available to the filter. One 
possible operation sequence that the controller may assume for a 
parallel arithmetic realization using an accumulator is listed below: 

1. Activate the A/D so that e^(kT) can be obtained. 

2. Shift the m(kT - T) and m(kT - 2T) storage registers to 
perform time delay. 

3. Clear the accumulator and its associated shift registers. 

4. Load ci and m(kT - T) into their respective shift registers. 
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5. Perform the multiplication c-^tnCkT - T) . 

6. Clear the accumulator shift registers. 

7. Load C 2 and m(kT - 2T) into the shift registers. 

8. Perform the multiplication C 2 in(kT - 2T). 

9. Clear the shift registers. 

10. Load ag and e^(kT) into the shift registers. 

11. Multiply age^CkT). The output of the accumulator now contains 

a^CkT) + c^m(kT - T) + c 2 m(kT = 2T) = eg(kT). 

12. Clear the accumulator and the shift registers. 

13. Load -b^ and m(kT - T) into the shift registers. 

14. Multiply -b^ m(kT - T) . 

15. Clear the shift registers. 

16. Load -b 2 and m(kT - 2T) into the shift registers. 

17. Multiply -b 2 m(kT - 2T) . 

18. Clear the shift registers. 

19. Load e^(kT) and 1.0 into the shift registers. 

20. Multiply l.Oe^(kT). The output of the accumulator is now 

e ± (kT) - b 1 m(kT - T) - b 2 m(kT - 2T) = m(kT). 


Repeat same process again. 

To ensure the correct sequence of operations that must be performed 
and the correct transfer of data such that there will be data available 
when required, the controller is divided into two parts: 1) the control 
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function generator and 2) the data transfer logic. The control function 
generator is simply a logic circuit which has as its output a sequence 
of pulses (with their timing and spacing very important) which controls 
the operation of the arithmetic unit, input/output, and memory. The 
data transfer logic is simply combinational logic circuitry which, under 
command from the control function generator, transfers data to and from 
the memory, input/output, and arithmetic unit. 

The first step in designing any controller is the selection of an 
operation sequence, such as the previous example. After this is 
complete, the control function generator and then the data transfer 
logic may be designed. 

Now that an approach has been presented by which all the functional 
components of a digital filter may be designed, the only remaining 
design consideration remaining is the interconnection of all four 
functional components such that they may function as a digital filter. 

The interconnection of the functional units may be done in numerous 
ways. No specific approach will be given here since, in most cases, 
each designer of a digital filter has what he thinks is his own unique 
and novel way to interconnect the functional components. In general, 
the connections of Fig. 4.8 must be made in as modular fashion as 
possible for possible LSI realizations. If they vary it will be because 
of programming form variation and order of D(z) variation. 
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Now that a basic SP computer implementation of a digital filter 
has been discussed, it will be in order to discuss variations of the 
SP computer implementation. 

Implementation by Microprogrammable SP Computer . 

The three previously discussed implementation techniques all 
have limitations, the most noticable of these being that each type 
implementation will realize a D(z) in only one programming form. It 
would be advantageous to design a digital filter which would realize 
a D(z) in any of the eleven previously discussed programming forms. 

In answer to the question that may arise as to why is one pro- 
gramming form better than the other; it can be shown, as discussed 
previously, that for a particular D(z) with set coefficients and 
sampling rate, different programming forms have different quantization 
errors. In general, for a particular D(z) that must be realized, it 
is desirable to choose the programming form that introduces the least 
quantization error, therefore necessitating the need for a digital 
filter implementation that can realize a D(z) in several programming 
forms. The implementation approach taken to do this is a microprogram- 
mable design as discussed in [12]. 

It will be the purpose of this section to give a discussion of 
the computer organization and not to go into too much detail about the 
logic design since it is not necessary for an understanding of the oper- 
ation of the microprogrammable implementation. 
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The prime function of the SP computer is the realization of a 
second-order digital filter in a choice of digital filter programming 
forms. As for the previous implementations this calls for the solution 
of the appropriate difference equations and entails the operations of 
addition, multiplication and time delay. The SP computer is binary, 
synchronous and parallel, with the calculations to be done using 2's 
complement arithmetic. It is a stored program computer, i.e., the 
Control Unit looks up the sequence of instructions in the memory, 
initiates arithmetic operations and causes operands or immediate results 
to be transferred between the Arithmetic Unit and the Memory, as 
required by the program instructions. Synchronization of the computer 
operations and generation of control pulses is the task of the Control 
Unit. 

To increase the sampling rate of the filter, a high speed Wallace 
Multiplier is employed as described in [12]* Also a rapid-access 
memory will be used for the program memory to aid in decreasing the 
worst case delay between the time of input sampling and its corresponding 
output. 

A block diagram illustrating the four basic functional units of 
a digital computer is as shown in Fig. 4.7. Therefore, just as for 
the previously described digital implementation of a digital filter, the 
organization of the SP computer will be divided into these four units. 

A detailed block diagram of the SP computer is shown in Fig. 4.19. 

The functional blocks in Fig. 4.19 which make up the Input/Output Unit 
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of Fig. 4. 7 are the A/D converter along with its input register and 
the D/A converter and its output register. In many cases the input 
register is considered as a part of the A/D converter and is not shown 
separately. The Memory is composed of the coefficient storage, variable 
storage, the program memory, and their appropriate memory address select 
and data select networks. The coefficient register, variable register, 
accumulator register, high-speed multiplier and the reduction logic 
network comprise the Arithmetic Unit. Included in the Control Unit 
are the sample clock, fundamental clock, binary counter, timing level 
generator, control matrix, the instruction register and its code 
translator. In addition to the four basic functional units which are 
generally used to represent digital computers, the SP computer in 
Fig. 4.19 employs an Operator Control Unit. The Operator Control 
Unit is used to program the SP computer for a particular filter form 
and also allows the operator to load the coefficients of the transfer 
function, D(z), into the coefficient storage locations. The Operator 
Control Unit can also be used in the testing and trouble shooting of 
the computer. 

At this point it would be desirable to present a brief description 
of the operation of the SP computer. For purposes of illustration, 
assume that a known transfer function, D(z) , is to be realized and that 
a specific filter programming form has been selected. With reference 
to the functional blocks of Fig. 4.19, the programming form is 
chosen by setting the 4 bits in the program register to the proper values 
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as will be defined later. Next the constants for the difference equations 
are manually written into the coefficient storage after switching the 
memory control switch, S, to the "1" position. Also, to control the 
output reduction logic, load the shift key words into the variable 
memory which are called by the quantize instruction to provide the 
desired number of accumulator bits to be input to the D/A or to be 
written into the variable storage. The programming of the SP computer 
is complete and it now awaits the first input signal. 

Normal operation begins as the sample clock gates the translated 
code of the program register into the program memory address register. 

The memory address register contains the address of the first instruction 
needed to calculate the difference equations for the filter programming 
form chosen. Stored in the program memory in groups of consecutive 
addresses are the macro -instruct ions for all filter programming forms. 

For example, the first seventeen locations in the program memory are 
the macro-instructions for the direct programming form. Upon receiving 
a read pulse, the program memory loads the 16-bit instruction register. 

The instruction format is shown in Fig. 4.20. Its four MSB's comprise 
the op code (operation code), the next six bits contain the first 
operand address (coefficient address), and the last six bits contain 
the second operand address (variable address). After the op code is 
decoded into one of the eleven available macro-operations, the Control 
Unit generates a corresponding sequence of micro-operations. Thus, 
each macro-operation is built up as a sequence of elementary micro- 


operations. 
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Op Code 


Second Operand Address 




16-Bit Instruction Register 


Fig. 4.20. Macro-instruction format. 
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The remainder of the discussion of a micro-programmable realizations 
will be devoted to the specification and organization of the four basic 
functional blocks of Fig. 4.7. 

1. Input /Output Unit . 

As before, it is the task of the input interface element to digitize 
the input analog signal e^t) in the A/D converter and furnish this 
digital signal e^kT) to the SP computer. The output, e 0 (kT), is in 
digital form and is converted to a discrete analog output, eg(t), by 
the output interface element, the D/A converter. 

A/D Conversion . The A/D converter will be such that it may trans- 
form e^t) into a two's complement binary representation and thus will 
be used as the input interface of the SP computer. Word size may vary 
from several bits up to sixteen bits including the sign bit. It is 
assumed that an A/D converter which meets the resolution and speed 
requirements of the SP computer is available. 

A 16-bit input register is used to hold the A/D converter output. 

The input register functions to maintain constant values for the input 
data bits during the period they are used. After each conversion the 
A/D converter produces an end-of-conversion (EOC) signal which is used 
to load the new digital word into the input register. If the A/D 
converter wordlength is smaller than sixteen bits, the remaining bits 
of the input register must be filled in with the sign bit or zeroes. 

D/A Conversion . Since the D/A converter is the output interface 
element, the selection of a D/A converter type is determined wholly by 
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the requirements placed on its output signal by the external system. 

This may lead to a D/A converter which converts a specified number of 
bits into either a unipolar or bipolar analog signal, the external 
system may require a pulse-width-modulated voltage signal, or even yet, 
may require a digital input, whereby a D/A converter is not needed. 

For the particular micro-programmable digital filter being discussed, 
the D/A was chosen so that a bipolar analog voltage may be presented 
at the output. 

2. Arithmetic Unit . 

It is the purpose of the Arithmetic Unit to perform the multiplica- 
tion and accumulation operations required to solve the difference 
equations of the digital filter programming forms described earlier. 

A major portion of the Arithmetic Unit consists of a high-speed 
multiplier. Inputs to the Arithmetic Unit are the multiplier (coefficient 
register), the multiplicand (variable register), and the accumulator 
register. Both multiplier and multiplicand inputs are 16-bits long 
and are coded using the two's complement representation. The output 
of the high-speed multiplier is a 34 bit, two's complement number which 
is stored in the accumulator register. 

High Speed Multiplier and Accumulator Register . The data input 
registers (coefficient and variable) and the data output register 
(accumulator) are organized to permit the multiplication of the contents 
of the input registers and to add the results to the contents of the 
accumulator. The simplified block diagram of the Arithmetic Unit in 
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Fig. 4.21 will be used to explain the calculation of an example 
difference equation. Note that the input to the accumulator register 
is gated. This is necessary since the output of this register is fed 
directly into the multiplier, i.e., the next state of the accumulator 
register is dependent not only on the multiplier inputs but also its 
own present state. Thus once the output of the multiplier reaches a 
stable state the data is gated into the accumulator register. 

The example difference equation requires two multiplications, 
with the results of each added to the contents of the accumulator 
register. Before starting the calculations, the accumulator register 
is cleared. Its initial contents are denoted by (Acc)q. Note that 
the accumulator register supplies an input to each column (which corres- 
ponds to a Wallace multiplier tree) of the array. This increases the 
size of the multiplier structure slightly but eliminates the time delay 
of an adder network which would other wise be necessary to add the 
contents of the accumulator and the multiplier output. After the 
first multiplication, the result is gated into the accumulator, be- 
coming (Acc)p For accumulation of the second product, a^e^(kT - T) , 
(Acc)^ and the partial products from the AND gate array are used as 
inputs to the free network. The result, 

e 0 (kT) = (12/16) (10/16) + (-10/16) (-6/16) = 180/256 

is gated into the accumulator register after sufficient time for the 
inputs to propogate through the Wallace multiplier trees. 
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Example difference equation: 


e 0 (kT) = a 0 e.(kT) + a^.OcT - T) 

let a Q = 0.75 = 12/16 = 0.1100 

a x - -0.625 = -10/16 = 1.0110 


e i 

(kT) - 0.625 = 10/16 = 0.1010 


e i 

(kT - T) = -0.375 = -6/16 = 1.1010 


0.1100 

1.0110 


0.1010 

1.1010 


0.00000000 

(Acc) o 0.01111000 

(Acc) ^ 

0.00000000 

0.00000000 


0.0001100 

1.1110110 


0.000000 

0.000000 


0.01100 

1.10110 


0.0000 

0.1001 


0.01111000 

(Acc) 1 = 120/256 1 



0.10110100 

(Acc) 2 = 180/256 


Fig. 4.21. 


Calcuation of example difference 
equation by Arithmetic Unit. 
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The first multiplication and accumulation operation is straight- 
forward since both the multiplier and multiplicand are positive numbers. 
However, this is not the case for the second operation. Note, that if 
the product of the multiplier bit and the multiplicand sign bit is a 
"one", it must be repeated in the left-most positions. In the first 
operations, this product was "zero" in all cases and thus the left- 
most fill-in positions contain all zeroes. Also, note that the 
negative multiplier in the second operation requires that the last 
partial product term be the two's complement of the multiplicand. 

This is accomplished by taking the one's complement of the multiplicand 
and forcing a "1" into the Wallace tree for the LSB of this partial 
product. This procedure is taken care of by the partial product 
generation logic of the high-speed multiplier and is presented in 
detail in [12]. 

Every phase of the multiplier has been demonstrated, except the 
establishment of the length of the accumulator register. Multiplying 
two 16-bit, sign two's complement numbers yeilds a 31-bit product. 

Since the direct digital filter programming form requires the summation 
of the greatest number of products (5) in the solution of any single 
difference equation, this sets the required length of the accumulator 
register at 34 bits. This means that if two maximum size 16-bit 
numbers are multiplied and accumulated five times, a 34-bit register 
would be required to express the sum. 
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Reduction Logic . There are several register lengths in the SP 
computer which need to be analyzed; the 34-bit accumulator register, 
the 16-bit data locations in the variable storage, and the N-bit D/A 
converter. This is brought about by the use of the solution of the 
difference equation in later calculations or as an output, e 0 (kT). 

In the first case, this 34-bit solution must be stored in a 16-bit 
location. Therefore it is necessary to quantize the output of the 
accumulator register to 16 bits. In the second case, it is also 
necessary to reduce the word length, since a 34-bit D/A would be both 
expensive and impractical. Thus the need for reduction logic has been 
clearly established. 

The number of bits in the output variable remaining after quantizing 
the contents of the accumulator register is determined by the size of 
the D/A converter. For an N-bit D/A converter the reduction logic may 
select the first N least significant bits (LSB's), the first N most 
significant bits (MSB's) or any intermediate group of N bits. 

, This selection of output data bits is accomplished by writing 
into the variable storage a shift key word for each quantize instruction 
in the program memory for a particular programming form. Part of the 
instructions will contain the address of this shift key, which is read 
from storage and gated into the reduction shift register, the contents 
of which determine those bits of the accumulator register to be loaded 
into the D/A or the variable storage. This is done by choosing the 
appropriate bits of the accumulator and shifting these bits to the 
left until they occupy those positions with output lines. 
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3. Memory . 

Each of the sections of memory shown in Fig. 4.19 provides 
storage for a specific type of data. The memory, as previously dis- 
cussed, is not only used for storage purposes, but is also utilized 
to perform the time delay operations which are required to calculate 
the difference equations. Included in the description of the memory 
sections will be the memory address registers, the address select 
logic and the data select logic. 

Coefficient Storage . Upon selection of a transfer function and 
a filter programming form, it is necessary to store the proper filter 
constants and coefficients to be used in the solution of the difference 
equations. The coefficient storage is a high-speed Read/Write memory 
which is composed of two memory modules. Each module has addressable 
storage locations for sixteen 8-bit words. Proper connection of the 
modules yields 9 (16 x 16)-bit memory as shown in Fig. 4.22. 

Data is manually written into the coefficient storage locations 
through the use of the data register, the memory control switch, manual 
address load pulse, and the manually controlled write pulse on the 
control panel. The address select logic in Fig. 4.22 is necessary 
to allow the coefficient address to be chosen from the control panel 
address register when manually writing coefficeint values into memory 
or the first operand address portion of the instruction register when 
reading coefficient values from memory. 
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Variable storage . The organization of the variable storage, as 
shown in Fig. 4.23, is similar to that of the coefficient storage. 

Its function is not only that of storing the input, e^(kT), and the 
output, e 0 (kT), but also that of performing time delay. Assume the 
input sample, e^(kT) and its previous value, e^(kT - T) , are stored in 
memory. After the last computation involving e^(kT - T) is completed, 
e^(kT) is written into the memory location allocated to ei(kT - T) 
where it waits until the next sampling period to be used as e^(kT - T) ; 
thus the time delay operation is performed. 

Note that the variable storage requires both data and address 
select logic. This is due to the variable storage being used by a 
multiple of sinks and sources. Data inputs to the variable storage 
may come from either the control panel, A/D, variable register or the 
reduction logic. The variable storage may be addressed by the control 
panel or the second operand address portion of the instruction register. 

Program memory . The program memory functions as a storage location 
for the macro-instruction necessary to solve the difference equations 
of the various filter programming forms. These instructions are 
organized into groups; each group corresponds to a particular 
programming form. Within each group, the macro-instructions are 
sequentially arranged as needed in the solution of the particular set 
of difference equations. 

The program memory is a high-speed ROM that has 256 words of 16 
bits each. Each word (macro-instruction) is broken into three sections: 
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the op code (4 bits) , the first operand address (6 bits) and the second 
operand address (6 bits). The op code specifies the operation to be 
performed on the operands in the memory located at the addresses 
specified by the two address fields. For op codes which require only 
one operand (or no operand) , that portion of the program memory is 
blank, i.e., contains any combination of zeroes and ones. An example 
is the store input instruction. Here the op code is "0101", the coef- 
ficient address is blank and the variable address contains a 6-bit code 
specifying the storage location in the variable storage which is to 
receive the input data. 

Since the program memory contains groups of macro-instructions 
for all filter programming forms, it is necessary to be able to locate 
the first address in each group once a form has been chosen. The 
contents of the program register is translated into the program memory 
address of the first macro-instruction for each filter programming 
form. Next, the sample clock gates the translated code into the program 
memory address register (PMAR) and the first macro-instruction is 
accessed by a memory read pulse. First and last in every sequence of 
macro-instructions is a start A/D instruction. This is necessitated 
by the overlapping of the instruction and execution cycle. 

In order to program the SP computer, it is necessary to code 
each digital filter form using 4 bits. Table 4.1 presents the coding 
scheme for the filter being described in [12]. 
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TABLE 4.1. 

Program Register Contents 
0000 
0001 
0010 
0011 
0100 
0101 
0110 
0111 
1000 
1001 
1010 


Program Code 

Digital Filter Programming Form 
Direct 

Modified Direct 
Standard 

Modified Standard 
Canonical 

Modified Canonical 

Parallel 

Cascade 

Modified Cascade 
Structure XI 
Structure XII 


1111 


Test Mode 
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As an example of the sequence of macro-instructions in the program 
memory, consider the direct programming form with the difference 
equation. 


e 0 (kT) = aQe^kT) + a^tkT - T) + a 2 e i (kT - 2T) 

-bjeofkT - T) - b 2 e 0 (kT - 2T) 


(4-18) 


Table 4.2 presents a word description of the macro-programming 
instructions for this form. 

Notice that during the time the A/D is converting the analog 
voltage signal to a digital signal the computer is calculating the 
portions of the difference equations that do not require the digitized 
input, e^(kT). This is intended to minimize the time delay between the 
input sample and its output response. 

The macro-instruction "Store (variable)" means load the contents 
(where "contents" is denoted by the enclosing parentheses) of the 
variable register into the specified address of the variable storage. 
Store (reduction logic) means to write the output of the reduction 
logic into the variable storage. In order to permit the use of A/D 
converters with different conversion rates, a "wait on A/D" instruction 
is incorporated. If the converison is complete at the time that this 
instruction is reached, the next instruction "Store input" is executed; 
otherwise, the computer idles until it receives the end of conversion 
signal form the A/D. However, there may be instances when a high-speed 
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TABLE 4.2. Program Memory Contents for Direct Form 


Op Code First Operand Address Second Operand Address 

Start A/D 
Clear Accumulator 


Multiply & Accumulate 

a 2 

e ± (kT - 2T) 

Multiply & Accumulate 

a l 

ei (kT - T) 

Store (variable) 


ei (kT - 2T) 

Multiply & Accumulate 

-b 2 

e 0 (kT - 2T) 

Multiply & Accumulate 

- b l 

eo(kT - T) 

Store (variable) 


e 0 (kT - 2T) 

Wait on A/D 



Store Input 


e i (kT - T) 


Multiply & Accumulate 
Quantize (Acc) 

Store (Reduction Logic) 
Quantize (Acc) 

Load D/A 

Reset RMAR & Halt 



Start A/D 
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A/D converter is employed. In this case the input e.j_(kT) is stored 
and not utilized in the calculations until all other multiplications 
are performed. This may also be an unnecessary delay in the response 
of the filter. Thus, it may be appropriate to postpone the "Start 
A/D" instruction in the program to ensure that the "Wait on A/D" 
instruction places the computer in an idle state. 

The codes for the macro-instructions are listed in Table 4.3. 

This table is used to generate Table 4.4, which presents the actual 
contents of the program memory, the variable storage and coefficient 
storage for the direct filter programming form. Addresses for the program 
memory are encoded in octal in Table 4.4. 

Note that in Table 4.4 the first and second operand address 
fields are 6-bits, while the actual memories (coefficient storage and 
variable storage) have only 4-bit addresses. Thus after loading 
the instruction register only the four right-most bits of each of the 
address fields are used as addresses for reading the contents of the 
coefficient and variable storages. Notice, also, the shift key word 
in the variable storage at location, 0100. This shift key is requested 
by the quantize op code, 0111, and is loaded into the reduction shift 
register. A shift key word exists for each quantize instruction. 

Space is allocated in the program memory for the remaining pro- 
gramming forms. They will not be illustrated since one may get an 
idea of their structure from observing the direct programming example. 
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TABLE A. 3. Macro-Instruction Op Codes 


Op Code 

Operation 

0000 

Start A/D 

0001 

Clear Accumulator 

0010 

Multiply & Accumulate 

0011 

Store (variable) 

0100 

Wait on A/D 

0101 

Store Input 

0110 

Load D/A 

0111 

Quantize (Accumulator) 

1000 

Store (Reduction Logic) 

1001 

Halt 

1010 

Reset PMAR and Halt 


Not Used 
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TABLE 4.4. Memory Contents for Direct Form. 


Program Memory 




First Operand 

Second Operand 

(Address) g 

Op Code 

Address 

Address 

000 

0000 



001 

0001 



002 

0010 

000000 

000000 

003 

0010 

000001 

000001 

004 

0011 


000000 

005 

0010 

000010 

000010 

006 

0010 

000011 

000011 

007 

0011 


000010 

010 

0100 



Oil 

0101 


000001 

012 

0010 

000100 

000001 

013 

0111 


000100 

014 

1000 


000011 

015 

0111 


000101 

016 

0110 



017 

1010 



020 

0000 




Address 

Coefficient Storage 

Address 

Variable Storage 

0000 

a 2 

0000 

ei(kT - 2T) 

0001 

a l 

0001 

ei(kT - T) 

0010 

-b 2 

0010 

eo(kT - 2T) 

0011 

- b l 

0011 

e 0 (kT - T) 

0100 

a 0 

0100 

Shift Key 

0101 

0101 

Shift Key 
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Program Memory Address Register . Macro-instructions retrieved 
from the program memory are read from memory addresses contained in 
the program memory address register (PMAR) . Inputs to PMAR came from 
the code translator and the control unit. Data from the code translator 
is loaded into PMAR by the sample clock, whereas, upon receiving the 
control pulse, "up date PMAR", PMAR performs the operation 

(PMAR) + 1 -*• PMAR, 

and now contains the address for the next macro-instruction. Since 
there are 240 addresses in the program memory, the PMAR must be eight 
bits long. 

4. Control Unit 

The program for a digital computer consists of a set of machine 
operations such as addition and multiplication. Instructions for these 
operations have been referred to as macro-operations. Inside the SP 
computer these operations are further decomposed into a set of elementary 
operations called micro-operations. Count, shift, gate the memory 
address register, are examples of micro-operations. Normally, in 
general purpose computers, macro-instructions are at the programmer's 
disposal and may be readily changed by re-writing the program. However, 
as previously described, in this SP computer the macro-instructions are 
pre-programmed in the program memory. These instructions are sequentially 
read from the program memory and loaded into the instruction register. 

In the op code portion of the instruction register is a 4-bit code 
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which specified the macro-instruction to be executed. This op code 
is fed to the control unit of the SP computer and a sequence of micro- 
operations is performed for each op code. 

It is the purpose of the Control Unit to translate the op codes 
and supply all synchronization and control pulses to the rest of the 
SP computer. A block diagram of the Control Unit elements is shown in 
Fig. 4.24. These elements include the instruction register, a decoder, 
a four-state counter, a timing level generator, a control matrix, a 
fundamental clock and an operation flip-flop, D. The organization and 
function of each of these elements is discussed below. 

Instruction register and decoder . As macro-instructions are re- 
trieved from the program memory, they are stored in the instruction 
register until the instruction is executed and a new one is retrieved. 

Only the op code portion of the instruction register is employed by the 
rest of the Control Unit, while the operand addresses are sent to the 
coefficient and variable memory address registers. 

There are eleven macro-instructions which are used by the SP computer 
to solve the difference equations of the various digital filter programming 
forms. Each op code portion of a macro-instruction is translated by 
the decoder into one of the macro-operations, f^, i = 0, 1, 2, *** 10. 

For each op code, the decoder activates one and only one output line. 

Fundamental clock . Employed as a fundamental clock is a free- 
running multivibrator whose frequency is obtained after determining a 
basic timing cycle for the SP computer, which can be done after the 



OP Code |First Operand | Second Operand 
~T Address Address 
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Fig. 4.24. Control unit. 



3-96 


detailed logic design of each functional unit is complete. Note that 
a clock signal that is 1/4 as fast as the fundamental clock is used as 
input to the Control Unit. This is because one of the micro-operations, 
quantize (Acc) , is a serial operation and requires these high frequency 
pulses in order to avoid slowing the operation of the SP computer. 

Four stage counter and timing level generator. There are fourteen 
micro-operations (from which eleven sequences of micro-operations are 
formed for the eleven macro-operations) , and the control signals for 
these micro-operations are designated as mj, j =0, 1, •••, 13. The 
number of timing levels for the execution of the macro-instructions 
will require six states of the four stage counter, L, for all instructions 
except multiply and accumulate, whose last micro-operation is performed 
on state fourteen. Thus this macro-instruction is used as a control 
input to the counter. 

Control matrix . Basically, the function of the control matrix 
is to provide the proper combination of op codes and timing levels. 
Essentially, .the control matrix is an AND-OR switching matrix. 

Operation flip-flop . Controlling the operation of the SP computer 
is the operation flip-flop, D. If D is set, clock pulses enter the 
four stage counter under normal operations. Flip-flop D may be set by 
three signals, the manual start, the end-of-conversion (EOC) signal 
from the A/D and the sample clock pulse. It may be reset by the control 
pulses, wait on A/D, halt, and reset PMAR and the manual stop button. 

Note that when the computer is halted by reset PMAR, the first instruction 
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to be executed is already in the N register ready for execution. Care 
must be taken to ensure that EOC does not occur before the wait on 
A/D signal during the execution of a given filter programming form. 

5. Operator Control Unit . 

As mentioned earlier, the Operator Control Unit is used in manually 
programming the SP computer to realize the transfer function, D(z), in 
the selected digital filter programming form. Referring to Fig. 4.19, 
the elements comprising the Operator Control Unit are the manual data 
register, the manual memory address register, two write pulse circuits 
(one each for the coefficient and variable storage), a program register, 
a memory control switch, a manual start pulse circuit, a manual stop 
pulse circuit, two pulse circuits for manually loading the coefficient 
and variable memory address registers, a row of sixteen indicator lights 
and a monitor switch, and a sample clock enable switch. 

Each of the registers may be implemented with single pole-double- 
throw (SPDT) switches for ease of setting and resetting. 

Data may be entered into either the coefficient or variable storage 
by setting the manual memory control switch to the "1" position, setting 
the switches of the manual data register to the binary data value, 
setting the address into the manual address register, loading this 
value into the coefficient or variable memory address register, and 
pressing the coefficient or variable write button. Notice the data 
register and memory address register are common to both the variable 
and coefficient storage since each employs a separate load address and 


write pulse. 
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Previous discussion has described the functions of the program 
register, manual start button and the manual stop button. The indicator 
lights and monitor switch function in checking the operational status 
of the SP computer. 

From the above discussion, one should understand the basics of a 
microprogrammable digital filter implementation. If further details 
are desired, one may refer to [12]. 

Next to be discussed will be the aspect of timesharing digital 
filter implementations . 

Time-Sharing of a Digital Filter Implementation 

For many applications of a digital filter it is sometimes necessary 
to provide discrete-time filtering for several independent loops or 
channels with examples being in control systems and digital communications 
systems. This may be accomplished with several digital filters but it 
would be more practical, more reliable, and more economical to provide 
this filtering with a special-purpose computer organized and programmed 
as a time-shared digital filter [9]. 

There are two basic ways in which a digital filter may be time- 
shared. The first of these is in multi-loops and the second is for 
higher order D(z) realization. 

Fig. 4.25 illustrates a digital filter being time-shared in N 
control loops. It is seen that there is an enormous amount of hardware 
saved by doing this. First there is only one computing element (SP 





Fig. 4.25. Block diagram of a digital filter being 
time-shared in multiple loops. 
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computer) required and second only one A/D is required, whereas if 
time-sharing weren't used N computing elements and N A/D converters 
would have been required. 

Fig. 4.26 illustrates a digital filter being time shared for 
a higher order D(z) realization. This can be done by cascade or 
parallel methods as shown in parts (a) and (b) of the figure. Cascading 
or paralleling filters to obtain a higher order realization is preferred 
because of the word length problem (resolution) that is encountered 
for a single high order D(z) realization. Time sharing of a single filter 
as shown in Fig. 4.26 can be easily accomplished if the assumption is 
made that n is even and that D(z) can be factored into n/2 second-order 
filters F^(z), ***, F ^( z ) an< * P artial fractional into n/2 second- 
order filters P 1 (z) , P 2 (z), ’ "^ n /2^ ' Each of the second-order filters 
may be realized by any of the available programming forms. The total 
transfer function D(z) still has the same memory requirements, the same 
number of memory locations in each filter storage unit, but now the 
variables and coefficients of each second— order filter are the signals 
stored. Since each second-order filter has the same program, the 
control sequence generator just repeats the same sequence n/2 times 
during each sampling interval. 

Both cascade and parallel techniques are desirable methods for 
realizing higher order D(z)'s, but the additional summing junction 
necessary with the parallel method makes the cascade method better 
suited for a modular realization. 
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For a better understanding of time-shared digital filter implemen- 
tations, let us look at the organization of a SP computer realization 
of a 2nd order time shared digital filter which may be time shared in 
two loops or cascaded or paralleled for the realization of a 4th order 
D(z) . 

Fig. 4.27 is a functional block diagram of the SP computer 
realization of two digital filters. Like previous SP computer realiza- 
tions, the computer functions are still divided into four main categories; 
input-output equipment, memory unit, arithemtic unit, and control unit. 

The input-output equipment provides the interface between the analog 
system and the digital computer. The analog inputs are sampled via the 
multiplexer and A/D, and the output samples are demultiplexed and shared 
in the buffered D/A's. 

The arithmetic unit, for this type realization, also must be able 
to perform multiplication, addition and subtraction. The selection 
of an arithmetic unit must first entail the selection of serial or 
parallel type arithmetic operation. Usually this choice means a decision 
must be made between computational speed and hardware economy. The choice 
made here was a reasonable fast and economical method where a parallel 
binary accumulator and two shift registers are used to add up the pro- 
duct terms as previously discussed. Serial arithmetic is used for 
partial product multiplication in the implementation discussed in [18] . 

As before, the reduction logic truncates to maintain word-length 
compatibility, and saturation logic sets the output word to maximum 
value when its range is exceeded. 
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The memory is divided into two identical storage units for filter— 1 
storage and filter— 2 storage. Corresponding filter signals such as 
e^(k - 1) (short notation for ei(k,T - T)) and e 2 (k ~ 1) have the same 
relative storage locations in their respective storage units and have 
the same addresses within their respective modules. The memory is 
controlled by Read and Write commands (not shown) and two address 
commands: a filter address to select the proper storage unit, and a 

signal address to select the proper signal location in a storage unit. 

The filter selection logic is controlled by the filter address and 
determines which storage unit is addressed by the signal address. This 
modular arrangement of the memory is .ideally suited for the "wired-OR" 
feature of integrated circuit memory. With this, feature storage for 
additional filters can be added by hard-wiring inputs and outputs 
of the storage units and simple modification of addressing. 

The control unit contains the master control unit that provides 
the multiplexing by means of the filter address and determines the 
sampling rate for each filter by controlling the start of each difference- 
equation computation. The control unit also contains the control-sequence 
generator which provides a sequence of instructions, initiated by a 
pulse from the master controller, to control the difference-equation 
computations. Thus the hard-wired program of the control-sequence 
generator determines the programming forms of the filters. To maintain 
simplicity, the same programming form is chosen for each filter, and 
hence the same sequence is generated each time regardless of which filter 




Fig. 4.27. Block diagram of the time-shared 

realization of two digital filters. 
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is being addressed and regardless of the number of filters being realized. 
A priority system is provided so that once the computation of a difference 
equation has begun, it is carried to completion even if it means that a 
sample for another filter must be omitted. 

Range Switching Digital Filter Implementation 

For many applications of digital filtering, it is desirable to have 
a very fast sampling rate. An example of this is when very high frequency 
signals are being filtered, and since the sample rate must be at least 
twice the highest signal frequency, it is seen that very high rates 
can be required. The limiting factors for a filter's sampling rate 
are the conversion speed of the A/D and D/A input-output system and 
the time required in the arithmetic calculations of the difference 
equations being realized. The most important factor which limits the 
speed of arithmetic calculations is the bit-lengths of the data processed 
internally by the filter. This includes the length of the input word to 
the filter, the internal variable wordlength (length of ej[(kT - nT) , 
m(kT - nT) , etc.) and the output wordlength. In general, the shorter 
these wordlengths are the faster the arithmetic calculations can be 
performed, with this being true for serial and parallel arithmetic 
type filters. Also, as a result of shorter wordlengths, the filter 
hardware is reduced. 

Because of the quantization error introduced by shorter wordlengths, 
a scheme must be devised to eliminate the larger quantization errors 
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Introduced by the shortened words. A method was devised to do this 
and the resulting filter was called a "Range-Switching" digital filter 
as described in [8], In short, the filter with reduced wordlength 
performed as if it had a much longer wordlength, under certain conditions. 

A block diagram of the filter described in [8 J is shown in Fig. 

4.28. Its operation will now be described. 

The A/D converter converts the analog input signal into an 8-bit 
digital approximation allowing the digital output of the converter to 
have an integer value ranging from 0 to 255. The input select logic 
selects either the four MSB's or the four LSB's of the A/D output to 
be the input to the SP computer. The scheme used is to select the four 
LSB's if all the four MSB's are logic "0" or the four MSB's if either 

one is a logic "1" for a signed magnitude code. It is seen from this 

that the input to the filter is in either of one or two ranges. If 

the input to the SP computer is from the lower range (four LSB's), 
it may have any integer value between 0 and 15 in steps of 1. If 
the input is in the higher range (four MSB's) it may have any integer 
value between 16 and 255 in steps of 16. Of the two ranges, it is 
seen that the higher range has the larger quantization step h and it 
is 16 times that of the lower quantization step. This also means 
that a four bit combination coming from the higher range input would 
represent a magnitude 16 times the same bit combination combing from 
the lower range input as shown in Fig. 4.29. Because of this, 
before the special purpose computer can process its four bit input each 
sample period, it must know the range from which it comes to properly reweight 
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Fig. 4.28. Block diagram of a "range-switching" 
digital filter. 
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the internal variables if a range change is seen from the last sample 
period. The internal variables of the SP computer are 8-bits long and 
can possibly be multiplied or divided by a factor of 16 each sample 
period before calculation of the difference equations start or their 
weight may remain the same. If there is a change in the input from the 
higher to lower range the internal variables are multiplied by a factor 
of 16 (the relative magnitude change of the input). It must be 
remembered that this is being done to keep the weight of the input 
relative to the internal variables. If the input represents a small 
value, the internal variables must be increased to make the input appear 
small. The opposite of this takes place when there is a change from 
the lower to the higher range. In this case the internal variables 
must be made to appear small to the large input, so they are divided 
by a factor of 16. If there is no range change between sample periods, 
the weight of the internal variables remains the same. The output select 
logic in Fig. 4.28 is employed to properly weight the 12-bit output 
of the filter. The output will have its greatest weight when the input 
for a particular sample period was in the higher range. If the input 
for a particular sample period is in the lower range, the output select 
logic weights the output bit configuration such that its analog voltage 
level representation is 1/16 of that of the same bit configuration with 
the input in the higher range. 

It was mentioned earlier that the "range-switching" scheme was 
used to reduce wordlengths but at the same time obtain the accuracy at 
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the filter output of a much longer wordlength filter. Also, the "range- 
switching" filter might be thought of as a technique by which for a 
fixed wordlength input, the quantization errors are reduced by the 
range-switching process . 

Filters of the "range-switching" design seem to have a bright 
future, especially with the advent of LSI. Using a reduced wordlength 
modular filter design, it would be very easy to have a digital filter 
composed of four or five LSI chips. 

One application of "range-switching" filters that has great promise 
is that of filtering in nulling type control loops. First, the design 
saves hardware because of the reduced wordlengths, it is fast and because 
of the design of the filter , -the quantization step length decreases as 
the loop error is nulled toward zero giving the loop finer granularity 
control to keep the error signal closer to zero. 

A block diagram illustrating the components of the SP computer 
realization of the "range-switching" digital filter is shown in Fig. 

4.30. 

The input/output unit consists of a successive approximation type 
A/D converter, a 12-bit D/A converter and the input/output select logic 
which are combinational circuits that select the input word and insert 
the output word into the proper bit position of the D/A. As previously 
mentioned the output of the filter is an 8— bit word. The 12-bit D/A 
converter is required such that the 8 output bits can be inserted into 
the 8 MSB positions of the D/A when the input is in the higher range 
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and in the lower 8-bit positions of the D/A when the input of the filter 
is in the lower range. In effect, a particular bit combination that 
is inserted into the lower 8-bit positions will have an analog voltage 
level 1/16 of what it would have been if inserted into the 8 MSB 
positions. It should be noted that the weighting factor of the output 
is the same as the input. 

Parallel fixed point arithmetic is used in the arithemtic unit 
in conjunction with shift registers, an accumulator, and the reduction 
logic to calculate the difference equations. 

The memory is composed of the coefficient storage (SPDT Switches), 
internal variable storage (flip-flop registers) , and the output storage 
(flip-flop register). The weighting of the internal variable memory is 
under command of the controller. Since fixed point arithmetic is used, 
multiplication or division by 16 is easily accomplished by shifting the 
internal variables left or right respectively four places relative to 
the binary point. If a weighted bit is lost when left shifting, provisions 
are made for the 8-bit internal variable to be saturated (all one's for 
a signed magnitude code) . 

In concluding the discussion of the hardware implementation, the 
control unit is composed of the control function generator and the data 
transfer logic. In addition to the normal tasks of the control function 
generator, it has the duty of deciding how the internal variables must 
be weighted each sample period. 
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A common question arises as to what is the limit on bit-length 
reduction. The most general answer to this is it depends on the appli- 
cation of the filter. The method of determining the minimum bit length 
is the trial and error technique of simulating the filter in whatever 
configuration it is to be used. As an example, the above described 
"range-switching" filter was used in a pendulous integrating gyroscopic 
accelerometer control loop. Time simulations in FORTRAN of the loop 
with the filter inserted demonstrated that the loop could be stabilized 
for pulse inputs to the loop with the filter having a 4 magnitude bit 
input 5 6 magnitude bit internal variables and a 8 magnitude bit output. 
From these bit lengths it is seen that an LSI realization of the filter 
would be quite small and simple. 

At the present time further work is being done to investigate the 
possibility of further bit length and hardware reductions of a digital 
filter such that LSI implementations will even be more attractive. 

LSI Circuit Digital Filter Implementation 

The first complete LSI implementation of a digital filter was 
designed and built by Autonetics Division of North American Rockwell, 
Anaheim, California. Concerning the state-of-the-art of digital filter 
implementation techniques, the design was the ultimate. The design 
is completely modular and therefore it is easily adaptable for an 
LSI realization. There are 2 main chips, a serial-parallel multiplier 
chip and a shift register chip<> 
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The serial-parallel multiplier chip is arranged such that it can 
perform all arithmetic functions required for a digital filter implemen- 
tation: addition, subtraction, and multiplication. The shift register 

chip contains shift register memories for the internal variables, input 
word, output word and also the control circuitry of the filter. 

The interface elements of the LSI implementation are external to 
the filter. The A/D output is fed into the filter in a serial manner 
and the filter output is also serial. The binary code used by the 
filter is two's complement. 

The internal wordlengths of the filter are adjustable. To change 
them, all one has to do is to make connection changes on the chips. 

The filter as designed has fixed length coefficients and each being 
16-bits long. They are easily set by single pole double throw switches. 

The filter realizes any first, second, or third order D(z) which 
has real poles in the parallel programming form. 

Commercial Digital Filters 

Now that several digital filter implementation techniques have 
been presented, a short discussion of digital filters available on the 
commercial market will be in order. Because of the relative newness 
of the area and the recent advent of LSI circuitry, there are few 
commercial builders around, two of which are listed below. 

One of the first builders of a digital filter for the commercial 
market was the Rockland Systems Corporation of Blauvelt, New York. They 
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now produce a line of programmable recursive digital filters which meet 
a wide variety of signal processing requirements. All of their filters 
are composed of four basic components - adders, multipliers, shift- 
register delays, and memory with a modular approach being adapted to 
provide the greatest possible flexibility and efficiency as is described 
in [18,19], The basic components are usually combined into second-order 
building blocks (two poles and/or two zeroes) and these blocks are then 
combined or multiplexed to realize any number of filters of any desired 
order. Programmability is achieved by employing a read/write coefficient 
memory. Fixed filter characteristics may be obtained with a read-only 
memory. Rockland produces a series of filters designed in the above 
manner. 

Rockland also produces a programmable tenth-order recursive digital 
filter which can realize arbitrary "all pole" designs such as Butterworth, 
Bessel, or Chebyshew low-pass, high-pass, or band-pass filters. Up to 
10 pole positions can be programmed through ten 12-bit filter coefficients; 
while up to 10 zero*s can be positioned at DC or the Nyquist frequency, 
or can be deleted altogether. Sampling rates of up to 100 KHz can be 
achieved. 

Electronic Communications, Inc. (ECI) is the second commercial 
builder of digital filters to be discussed [26]. 

The filters produced by ECI are the nonrecursive type and they are 
actually signal-processing instruments that perform the convolution 
integral in order to effect a filtering function. This means the filters 
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work strictly in the time domain and is represented by its impulse 
response and not by its amplitude and phase characteristics as are some 
nonrecursive filters. The analog input signal to this filter must be band- 
limited, and is accomplished by using a low-pass analog prefilter. The 
input bandwidth is restricted to less than half the sampling rate. One 
of their more common filters with a sample rate of 10 KHz has a pre- 
amplifier cut-off at 2.5 KHz. 

The ECI digital filter has two identical shift register memories. 

One stores samples of the band-limited input signal, while the other 
contains samples of the impulse function of the filter that is desired. 

The sampled impulse response is represented by coefficients that are 
the various amplitudes of samples spaced equidistant along the time axis. 

The coefficients are obtained by using computer software available 
from the company and programmed into the digital filter via a paper 
tape. The tapes are set up so that the programs can be put on a time- 
shared computer system. Up to 200 coefficients can be stored in the 
sampled-impulse-response memory. 

The contents of the two memories are fed into a single multiplier 
section that forms the product of corresponding samples form each 
memory. The output of the multiplier goes to an accumulator that puts 
out the digitized filtered version of the input waveform. 

The digital filter designed in this manner can only simulate zeroes 
for the filter transfer function because of its nonrecursive nature. This 
does not seem to limit its use though, as any filter can be represented 
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by its impulse response. There are few recursive filters whose impulse 
response cannot be satisfactorily represented by ECl's nonrecursive filter. 

This concludes the discussion on SP computer implementations. FFT 
hardware realizations will now be discussed. 



V. FFT HARDWARE 


We have previously seen that a digital transfer function, D(z), 
may be calculated by the FFT. The theory behind this was discussed, 
enabling the discussion of FFT hardware to now follow. Three area’s 
of FFT hardware will be discussed starting with commercial FFT pro- 
cessors that are available and concluding with a discussion of the 
MIT fast digital processor. 

Commercial Equipment 

There are several manufacturers of FFT processors. One of these 
is the Raytheon Computer Co. of Santa Ana, California. This company 
manufactures what it calls an "Array Transform Processor." This 
processor has several capabilities, among them FFT and Inverse FFT 
processing, convolution integral processing, complex multiplying, 
complex spectral magnitude processing, real multiplying, read add/ 
subtract processing, scanning of arrays and array movement. 

The Raytheon array transform processor is an auxiliary processor 
of the Raythean 700-Series computers and is a hardware array processor. 

It comes in several models with the models differing in the number of 
data points that can be handled (256 -*■ 16384) and the wordlengths available. 

Another company which manufactures a hardware FFT processor is 
the Elsyter Co. of Syosset, New York. Their processor is labeled as 
the 306/HFFT. It will calculate the direct or inverse FFT and it is 
contained within the mainframe of its host NOVA 800 computer for more 
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efficiency and small size* It has a core memory of 4096-16 bit words 
expandable to 32768 words* It has the features of hardware multiply 
and divide, teletype interface, array complex coordinate converter to 
perform Cartesian to polar and vice-versa conversions without program 
intervention and an hardware FFT interface subroutine to control the 
operation of the hardware FFT. 

It can be used in three modes: 1) Stand alone peripheral FFT 

processor, 2) part of a free standing spectrum analyzer system, and 
as 3) a free standing computer. 

The last commercial FFT processor to be discussed will be the "FFT 
256 FAST FOURIER TRANSFORM ANALYZER," a "stand-alone" processor manu- 
factured by Unigon Industries, Inc., Plainview, Long Island, New York. 
It is a smaller processor than the one previously discussed in that it 
has a capacity of 256, 1024, or 4096 real points and a wordlength of 
8-bits. 

This processor performs the functions of the direct and inverse 
FFT, power spectrum and cross power spectrum analysis, square spectrum 
analysis, auto correlation, cross correlation, convolution, convolution 
spectrum and auto correlation. Each one of these functions is switch 
selectable. Let us now discuss a fast digital processor designed 
by MIT. 

MIT Fast Digital Processor [27 ] 

There are many techniques by which GP digital computers may be 
modified or enlarged to increase the operating speed. As an example, 
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speed savings may be attained by attaching fast multiply and divide 
hardware or by having separately addressable memory modules so that 
instruction cycles and data cycles may be overlapped. Increases in 
speed results from attaching arithmetic hardware which performs high 
speed special operations such as digital filtering and discrete spectrum 
analysis. If this arithmetic hardware is added in addition to a high 
speed memory, speed increases on the order of 40 to 100 can be attained 
in performing operation such as the FFT. This technique has a disadvantage 
in that it requires programming that is not easily structured, which 
in turn decreases the speed advantage of the special hardware. 

It was because of this that emphasis was directed toward incorporating 
more general purpose features into a signal processing computer structure. 
What resulted was the MIT fast digital processor which is a general 
purpose digital attachment to a UNIVA/C 1219 computer. The architectual 
changes required to increase the speed of repetitive arithmetic operations 
for signal processing can be classified as 1) the use of scratch pad 
memories, 2) pipeline schemes, and 3) parallel processing. 

The fast digital processor (FDP) , designed with the above architectural 
changes in mind, is able to perform signal processing simulations close 
to two orders of magnitudes faster than present conventional digital 
computers. As an example, a vocoder simulation which normally requires 
about 200 times real time on a standard computer, could be programmed 
to operate close to real time on the FDP. 
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The main applications of the FDP are in the areas of communication, 
radar, speech processing, biology, medicine, and sonar. 

Fast commercial integrated circuit elements and a logical structure 
which permits each main unit of the machine to operate at maximum 
speed enables the FDP to obtain a speed advantage. 

The arithmetic section is designed to perform efficiently the 
sum-of-products operations which are most important to recursive 
digital filtering, the FFT and correlation operations. 

The data memories are structured so as to exchange data with the 
arithmetic section at maximum efficiency. 

The control uses a separate memory for storing instructions. Its 
structure allows the data memories to operate at maximum speed. 

Let us briefly describe some of the main features of the FDP 
structure. 

FDP structure . Fig. 5.1 shows the most important data transfer 
paths of the FDP. Programs are run from the memory M c , which controls 
the main data flow from memories M 3 and M b to all elements shown in 
Fig. 5,1, Since M c is intended to be a small memory, longer programs 
can be stored in M a and M b and block transferred to M c when needed. 

M 3 , M b , and M c are addressed independently and can therefore be operated 
in parallel. Except for block transfers from M a and M^, the control 
memory M c cannot be written into making the programs run by the FDP 
almost non-self-modifying. 
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The parallelism inherent in the FDP is partially indicated in 
Fig. 5.1. Listed below are a few special features incorporated for 
speed: 

1) four arithmetic units, each including a'multiplier which can 
operate in parallel with the main arithmetic registers; 

2) two independently addressable integrated circuit memories, 

M a and M* 3 with read and write times of 140 ns; 

3) a separate instruction memory, which allows overlap of instructions 
and data cycles; 

4) a double length instruction word which enables two instructions 
to be simultaneously executed on the FDP. 

The size of M a and is 1024 words and is addressed by a 10-bit 
word. Addressing is indirect, through M* 1 , a 16-register 24-bit integrated- 
circuit memory, as is shown in Fig. 5.2. The indexed address for M a 
and is formed by adding the contents of the two 12-bit portions of 
to X a and X* 3 . Writing into requires no special instructions 
because addresses 0 through 15 of M a and are wired to control 
as well as the data memories on a write cycle. 

The I/O capabilities of the FDP were minimized deliberately since 
the UNIVAC 1219 already supplied most of the necessary I/O control. An 
A/D and D/A converter were applied to make the FDP applicable for real- 
time processes. The only other I/O path is to and from the mother 
computer, the UNIVAC 1219, which will be needed for assembly ing and 
editing FDP programs and supplying medium size core storage. 
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The FDP is an 18-bit fixed point processor. Floating point routines 
may be used but must be programmed, as must multiple processor arithmetic. 
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